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BCA and HRA: Two Programs for 


Computing Economic Equilibria 


by 
Thomas Elken 


I. Introduction 


BCA and HRA are computer programs designed to solve a version 
of the economic equilibrium problem. For an introduction to this 
type of problem and models for dealing with it, the reader is directed 
to Elken [2]. Very briefly, we will describe mathematically the 
problem which these codes solve. 


Find x, t, A, and Cf such that 


oO, ie#l1, ..., 


Oo, i= IH +1, cers NROW 


where D is a matrix with NROW rows and NCOL columns, IH is 
the number of consumers or households, eNcoL is the NCOLth unit 
vector, and (+, +) is the usual inner product. 

The BCA is an implementation of the bilinear complementarity 
algorithm presented in Wilson [7] and in Elken [2]. The HRA code 
implements the homotopy retraction algorithm as described in [2]. 
The reader must consult the latter work to be able to formulate a 
problem correctly and to be cognizant of the conditions under which 
these algorithms will solve his/her problem. 

The algorithms have so many features in common that many of 
the subroutines in the two programs are identical. In particular, 


both algorithms begin by solving the linear program 


maximize XNROW 


subject to [IID] (=) =b, (2) 
t, 220. 


To solve this linear program, the code LPMl written by John 
Tomlin is used. Actually, the subroutines comprising LPM] are part 
of both codes, BCA and HRA. LPMI has been documented in J. Tomlin 
[4]. LPMl1 has been converted into a linear complementarity code, 


LCPL, which has been documented in J. Tomlin [5] and [6]. Since 


LCPL uses many of the subroutines of LPMl without change, the reader 
is referred to the works cited for a description of how the matrix 
A= [I|D] is stored and how transformations are processed. 

Next we discuss some important details of the implementation 


of BCA and HRA. 


Source Language 
The programs are written entirely in FORTRAN IV for the IBM 


series 360 and 370 computers. They are WATFIV compatible. 


Specification 


The main program is executed as a job step. The program 
input is described in Section II.1. All input occurs in subroutines 


MAIN, INPUT, and BASINP. 


Error Indicators 
Error indicators and other diagnostics and messages will be 


indicated in the documentation for each subroutine. The error indi- 


cators for the LPMl portion of the code are described in [4]. 


Subroutines 


The routines making up BCA and HRA are as follows: 


wceeie 


A OSA RCA A ARPT 


i BCA HRA 


MAIN MAIN 
BLCONS BLCONS 
PEND ics 06 pepe Seal eset) enafie ohare ta, “SMP 
BEVOE aes eer esl es aol ie a CR EVOD 
UPARO © oie twee ie 6 ss st ye, es 6 RARE 
BSONG. soi 6s ge Sow we ee ee ee BSENG 
SUPERB. <5 6 Woe ee ew we UPR 
RECALC . . 2 ss © © © © © © © © © heh e)6URECALC 
ENDPNT ENDPNT 
QUADS ° See. wee QUADS 
DSENT . . . DSENT 
GTN 
DERIVG 
CONCHK 
FIN 
DERIV 
NORM 
DECOMP 
SOLVE 
SING 
DEBUG 


wey ee 


The dotted lines indicate subroutines which are essentially 


identical. 


The remainder of the subroutines are from LPMl for both programs. 


ee 


BLOCK DATA 

INPUT 

FTRAN CHSOL 

BTRAN UPBETA 

FORMC NORMAL 
[ TRUEDJ ITEROP 

PRICE UNRAVL 

CHUZR CRASH 

WRETA BASINP 

SHIFTR CLEAR 

INVERT 

UNPACK 

SHFTE 


These subroutines are described in greater detail in the 


remaining sections of this report. 


Program Size 


The total length of BCA is 3831 source statements (including 
comments and common blocks, etc.); HRA is 3772 statements long. 


LPM1 comprises 1819 statements of each of these. 


Array Size 


The standard version of the codes requires a blank common 
array of 133,256 bytes of core. This is dependent on the setting of 


f maximum problem dimensions. 


There are eleven labelled common blocks which require 97,448 
bytes of core. Only the size of common blocks LNCONS and INDXZ 


are dependent upon problem dimensions. 


Ac curacy and Conver gence 


The linear constraints (1) and (2) are satisfied with a toler- 
ance of ZTOLZE with a default of 1.E- 4. The bilinear constraints 
(3) are satisfied to a tolerance of TOLFZ; 1.E - 10 is the default. 

Convergence is guaranteed in the BCA and HRA when the first 
IH values of b are chosen appropriately. See Elken [2] Chapters 
V-VI for a discussion of the convergence properties of these algorithms 


and suggestions for improvements. 


Timing 
Using the FORTRAN H compiler with OPT = 2. the codes have 


solved a problem with 6 consumers and 4 goods with D€ eel in 


-45 and .38 seconds for the BCA and HRA, respectively. A problem 


with 3 consumers and 56 goods and D €& pox 107 


required 4.54 and 6.19 
seconds for the BCA and HRA to solve. A description of numerical 


results for these and several more test problems is in Chapter V of 


tZ}. 


— 
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II. The Bilinear Complementarity Algorithm (BCA) 
II.1. Input Requirements 


The program input consists of 4 (or 5) segments. 


A) Parameter specification relevant to equilibrium problem 

B) Parameter specification relevant to the auxiliary linear program 
C) Data for the LP in MPS standard form 

D) Basis input (optional, see [5]) 

E) Specification of the C matrix (3). (Initial endowments for 


the individual consumers.) 


We will now deal with these in succession. 


A. Program Parameter Input 


Parameters are changed from their default value for each problem 
by means of a FORTRAN "NAMELIST" input from the card reader. The 
first card must start with &PARMl in column 2 or later; each card 


must begin with a blank, and the list must be terminated by &END. 


; Example 
columns 
£23 
&PARM1 IH = 3, L = 2, ICECHO = 1, &END 
The parameters which can be defined by the &PARM1 statement 


are listed below along with their default values. 


alee ee accede 


Cad 


A EOE, SNL 


REAL *4 Tolerances 
TOLFZ (DEFAULT = 1.E - 10) 
Termination criterion for the solution of the bilinear equations. 
TOLBD (DEFAULT = 1.E - 4) 
Feasibility tolerance for the linear constraints. 
TOLCV (DEFAULT = 1.E - 5) 


Termination criterion for the return to the curve. 


INTEGER *4 Parameters 

IH Number of consumers (traders). This parameter has no default. 
It must be specified in the &PARMl1 statement. 

ICNTRL (DEFAULT = 1) 
If ICNTRL-LT.O, then the program behaves just like LPMl1 to 
solve one or more linear programs. This may be useful when 
one desires to solve some preliminary linear programs to deter- 
mine initial utility levels. Otherwise, an equilibrium problem 
will be processed by BCA. 

IECHO (DEFAULT = 0) 
If IECHO = 1 while ICNTRL.LT.O, the D matrix of (1) 
will be printed out, ten columns at a time. 

ICECHO (DEFAULT = 0) 
If .FQ.1, the C-matrix of (3) will be printed. 

IPROD (DEFAULT = 0) 
If .EQ.1, it is assumed that the consumers own certain shares 
of the producing firms. Thus, data containing those shares 


must appear following the LP data in 10F7.4 format. 


ITLIME (DEFAULT = 500) 
The maximum number of cells which we allow the program to pass 
through. 
(DEFAULT = 0) 
Unit number of output of the optimal basis for the auxiliary 


linear program. If KOUTB.EQ.0O, no basis will be saved. 


REAL ¥4 Parameters 
STPMX (DEFAULT = 100) 

The maximum possible step in the direction tangent to the curve. 
STPRD (DEFAULT = 0.5) 

The factor between O and 1 used to reduce the steplength 


when it exceeds STPMX. 


B. Linear Program Parameters 

Next, a list containing parameters of importance for the linear 
programming portion of the code must be included. The beginning of 
this list is &PARAM. We leave the description of these paramaters 
to [4] with one exception. The parameter I0OBJ contains the row 
number of the objective row in the D matrix. Since this row is 


always (2) no flexibility is necessary here so we require that 


Pi 
NROW 
IOBJ be the last row of D. Since we did not want to alter the LPM1 


part of the code, the user must know how many rows D has and enter 
the appropriate value for IOBJ in the parameter list. This rather 
irritating requirement makes bookkeeping easier in many parts of the 


code. 


Also, the auxiliary linear program (4) is a maximization problem 
while LPM1 is a minimization code by default. This can be remedied 


in one of three ways: 


BY 
1) Let the last row of D be ~eyRow’? 
2) Set ZSCALE=-1.0 in the &PARAM list, or 


3) Set IOBJ = -NROW in the &PARAM list. 


Remember to do only one of the alternatives above. 


C. Problem Input 


The D and b_ referred to in (1) are read from unit KINP 
(DEFAULT = 5) in slightly modified "MPS format". The card images 


required are: 


NAME Card 
This has "NAME" in columns 1 to 4 and the (up to 8 characters) 
problem name in columns 15-22. This card is optional but highly 


desirable. 


ROWS Card 


Has "ROWS" in columns 1-4. 


Row Names 
Each row is assigned a (unique) name of up to 8 characters, 


one per card, in columns 5-12. Embedded blanks are allowed (unlike 


— 


MPS). Also row types must be supplied in columns 3 of each card 


L: < constraint 


o 
Vv 


constraint 


peo] 
i 


constraint 


N: objective function row . 


For the equilibrium problem all constraints will be of type 


"L" except for the objective function. 


COLUMNS Card 


Has “COLUMNS” in columns 1-7. 


Matrix-Element 


The non-zero elements of D are supplied by column. All 


the elements of a column must be together. Each column is assigned 


a (unique) name of up to eight characters. The format is: 


cols 1-4 5-12 15-22 25-36 40-47 50-61 
blank column row element second second 
name name value row element 
name value 


(optional) (optional) 


Again embedded blanks are allowed in column names. 


li 


RHS Card 


Has "RHS" in columns 1-3. 


Right Hand Side Elements 


The right hand side vector (b) may be given a name, which 
should be different from any row or column name. The elements are 


given in the same format as the matrix elements. 


ENDATA Card 


Has "ENDATA" in columns 1-6. Sample input is shown in Section 


D. LP Basis Input 


If KINB is initialized to a nonzero value in &PARAM, a basis 
is read from unit KINB. If KINB.EQ.KINP, the basis must be entered 


following the LP problem input, and before the ENDATA card. This basis 


will be the initial basis for solving the auxiliary LP (2). 


NAME Card 
Has "NAME" in columns 1-4 and the basis name in columns 
15-22. If this name does not agree with the problem name, or the 


card is missing, a warning is given. 


Basis Cards 
For basic columns of D, the column names are given (one per 
card) in columns 5-12. For basic columns of I (slack variables), 


the corresponding row name is given in columns 5-12. 


ENDATA Card 


Has ENDATA in columns 1-6. 


E. C Matrix Input 


In most problems the data for C will be predominately extracted 
from the RHS vector b. For this reason we use a rather short, but 
complicated, form for the input of C. 

The data is read from unit KINP directly following the MPS 
format data. The instructions constructed from the pointers in this 
data produce a matrix packed into a vector in the same manner that 
{I|D] is packed (see (6]). 

In the current version of the code it is assumed that (if 
IPROD.NE.O) each consumer I is associated with a fraction SHR(I), 
which represents his share of each firm in the economy, such that 
ash SHR(I) = 1.0. If this assumption is not true, the correct C 
matrix can still be input; however, it requires more work. 


The form of the data for inputting the C matrix is as follows 


1. One or more cards containing (SHR(I), I = 1, IH) in 10F7.4 


format (only if IPROD.NE.O). 


2. Several cards describing the column of C associated with 
consumer 1. 
3. A blank card. 


4. Information associated with consumer 2, etc. 


We now elaborate on the specific requirements of point 2 
above. Assume that we are dealing with column I of C. 

The first card for column I contains values for the integer 
variables ISHR, LL, and KK in 314 format. ISHR is a control 


parameter which behaves as follows: 


A) If ISHR.EQ.1, the program reads one or more consecutive ele- 
ments of the RHS corresponding to constraints on one or more of 
the firms in the economy. These numbers are multiplied by SHR(I), 
the Ith consumers share of the firm, and stored in C. 

B) If ISHR.EQ.0, the program reads one or more consecutive elements 
of the RHS corresponding to the constraints on consumer I alone. 

C) If ISHR.EQ.-1, the program reads one or more elements from the 


card(s) following in the input stream in 10F7.4 format. These 


usually correspond to commodities (or firms) which are owned by a 


number of consumers, but not in any fixed proportion. 


LL and KK determine the range of rows of C which will be 
altered due to this card. If only one row is involved, either KK 
can be set equal to zero (leave columns 9-12 blank), or set KK equal 


to LL. 


1 | 


A blank card indicates that the program should move on to the 
next column of C (the next consumer). 

If the problem has been formulated and input correctly, the 
sum of the columns of C must be equal to the RHS vector (see [2], 
Section V.3). Examples of input and output for some small problems 


will be given after the program in Section II.4. 


I1.2. Main Program 


The main program consists of 350 statements which perform a 
variety of functions: the parameters are initialized and altered 
by reading the name list PARMI, the subroutines of LPMl are called 
to solve the auxiliary linear program, the data is read so that the 
© matrix can be constructed, the machinery is set up for the bilinear 
complementarity algorithm to begin, the main path-following subroutine 
ENDPNT is called, the updates and decisions made necessary by the 
output from ENDPNT are performed, and the final solution is printed. 
If options were added so that a number of equilibrium problems could 
be solved, this main program would undoubtedly become sort of a master 
subroutine, as NORMAL is for LPM. 


Restrictions relevant to the use of BCA. 


1. The number of consumers (IH) must be less than or equal to 10. 


2. A must have not more than 350 rows or 400 columns. 


Of course these limits can be extended by redimensioning all 


vectors of size 10, 350 or 400 and allocating more core. 
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Before describing the subroutines of BCA we will define the 
variables in the labelled common blocks and some of those in blank 


common (those of use in the subroutines outside of LPM1). 


Variable Glossary 


We have specified that almost all real variables are REAL*8 


by the statement 


IMPLICIT REAL*8 (A-H, 0-Z) 


Notable exceptions are the initial data in A and C which are stored 


as REAL*4 variables. 


COMMON/LP1/: 

PI(1302) The complete set of dual variables (prices and relative 
costs). 

XX(1302) The complete set of primal variables, slacks first and 


then the structural variables in their original order. 


COMMON/BLCST/ (Bilinear constraints) 


COMMON/BLCST2/ 


These variables store the current version of the budget surplus 


function in terms of the superbasic variables (see Elken [2] Section 


III.3). The functional form is 


” 


Type 1: f. = BF1 + D1*PI(MU) - DIAG(PI(MU))*(F1*XX(NU) + El) 


(3) 


Type 2: f= BF2 + D2**PI(MU) - DIAG(XX(NU)) *(F2*PI (MU) + E2) , 


where PI(MU) refers to those components PI(J) of PI such that 
MU(I) = J. for some I = 1, 2, ..., KMU, KMU is the current number 


of dual superbasic variables, and 


PI(MU(1)) QO «ee 0 
0 PI(MU(2)) 0 
DIAG(PI(MU)) = : 
0 ° PI (MU (KMU) 


XX(NU) and KMU are defined similarly. 
Also in COMMON/BLCST/ are C(800), IC (800), and LC(20) 


which store C in packed form in exactly the same manner as A. 


COMMON/LNCONS / 


These matrices and vectors allow us to express the basic vari-~ 
ables in terms of the superbasic variables. If IH and DIGMA are 
index sets for the primal and dual variables, respectively, then the 


relationship is 


ALTER 


KNU 
XX(JH(I)) = J} GL(1,J)*xXX(NU(J) + BACT) 
J=1 


KMU 
PI(DIGMA(I)) = J} G2(1,J) *PI(MU(J)) + BB(1) 2 
J=1 


/ where the arrays are dimensioned 61(350,10), G2(400,10), BA(350), 


BB(400) . 


COMMON/INDX1/NUH(10), MUH(10), NU(10), MU(10) 


NUH(L) 


MUH (1) 


NU(I) 


MU (1) 


These arrays correcpond to 0, o, Y and Y in the description 


O if XX(I) is not superbasic 


if XX(I) is the Jth 
superbasic variable I= 1, .-.., IH 


if PI(I) is not superbasic 


is PI(I) is the Jth 


dual superbasic variable, I = 1, ..., TH 


if the Ith superbasic variable is XX(J) 


if the Ith dual superbasic variable is PI(J). 


Note that NUH(NU(I)) = I. 


COMMON/ INDX2/ JH ,DIGMA ,KINBAS , IDBAS 


of the algorithm in Chapter III of [2]. 
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JH(1) = J if XX(J) is basic and pivots on row I, 
DIGMA(I) = J if PI(J) is dual basic and corresponds to the Ith 
row of G2, 
O if XX(I) is non-basic 
KINBAS (I) = 
J if XX(1I) is basic and pivots on row J. 
Q if PI(I) is not dual basic 
IDBAS(I) = 
j if PI(1) is dual basic and corresponds to the Jth 
row of G2. 
COMMON/SCAL/ 
BT = A scalar that helps to define the normal hyperplane subproblem 
in ENDPNT. 
JI = The index of the constraint in Gl or G2 which defines the 
exiting facet for the current celi. 
MFLAG =A flag which defines the type of jacobian which must be 
calculated. 
DD1 = The number of bilinear constraints which are currently at 
zero. 
P = The index of the pivot column determined by subroutine FINDP. 
PD = The flag which determines whether the initial facet corresponds 
to a primal (PD = 1), a dual (PD = -1) variable, or a bilinear 
constraint (PD = 0) at zero. 
MPD = Identical to PD except that if refers to the final facet of 


the currant cell. 
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ain ibe tr 


INEQ 


KFUN 


KJAC 


COMMON/DIM/ (DIMENSIONS) 


IH 


M 


N 


DD 


BLL 


BL2 


1, if the bilinear inequality is the last of the Type 1 


constraints. 


2, if the bilinear inequality is the last of the Type 2 


constraints see (1). 


Counts the number of functionals evaluated in (1). 


The number of partial derivatives of the functions defined 


by (1) which are calculated during the course of the algorithm. 


= The number of households. 


Number of rows in the problem. 


Number of structural variables in the problem. 


Current number of superbasic variables (= DD1 + 1). 
= Current number of primal superbasics. 
= Current number of dual superbasics (DD = KMU + KNU) 


KMU - 1, if INEQ = 1 


KMU, otherwise 


KNU - 1, if INEQ = 2 


KNU, otherwise 


ae i : 


N+M 


en ARERR SN 


—_— = Ge 


—d 


COMMON/INT/ 


IPS (30) A vector storing the permutations for a LU decomposition 


linear equation solver used in conjunction with Newton's 


method. 
KDET The sign of the determinant of the current jacobian of f. 
KOUNT The number of times ENDPNT has been called. 
ISING A flag which is 1 when the current jacobian is singular, 


QO otherwise. 


COMMON/TOLER/TOLFZ, TOLBD, TOLCV, ICNTRL, IECHO 


EES3s 


(2) 


These parameters are defined above. 


Subroutines of BCA 
SUBROUTINE BLCONS: Calculates the coefficients of the first DD 
bilinear constraints. DDl of these are constraints which define 
the curve. The last is an inequality which is not yet binding, 
but which must be the next one to become binding. 

The calculations are motivated by the fact that the bilinear 


functions 


M 
} C(I,J)*PI(J) - PI(I)*XX(I), ral, asx DO 
j=1 


can be reduced to functions of the superbasic variables by using 
(6). The position of DD in the index sets MU and NU deter- 


mines BL1, BL2, INP, INQ, and INEQ. 
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(3) 


(4) 


(5) 


(6) 


(7) 


SUBROUTINE FINDP (PD1, IS, P): Finds the element of largest 
absolute value in the row determined by (PD1, IS) and stores 
the index of the variable corresponding to that element in P. 
When we.say the row determined by (PDI, IS) we mean GI1(IS,-) 
if PDl=1, and G2(IS,-) if PDI =-1. 

SUBROUTINE PIVOT (SS,RR): If variable XX(SS) is entering 

the basis and variable XX(RR) is leaving the basis, the columns 
of Gl and G2 andthe columms BA and BB _ must be updated. 
This is accomplished by a simple pivot. The details for this 
pivot are contained in Section III.5 of [2]. 

SUBROUTINE UNPACKC(II): A call to this subroutine causes the 
IIth column of C to be unpacked and stored in the M-vector Y. 
SUBROUTINE BSCNG(S,R): This subroutine updates the index sets 
JH, KINBAS, DIGMA, and IDBAS when variable S replaces variable 
R in the primal basis. 

SUBROUTINE SUPERB(KEY, PD1, IS, PD2,JS): Revises the index sets 
defining the superbasic variables. It also adds or removes 


columns of Gl and G2 depending on the values of the parameters 


0, if columns are both added and removed, 


KEY = {1, if columns are only added, 


2, if columns are only removed. 


(PD1,1IS) determine the column to be added. 


(PD2,JS) determine the column to be removed. 
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The calculations iuvolved in calculating the column to be added 
are described in [2] Section III.3. 
SUBROUTINE RECALC: The parameter INVFRQ of LPM1 specifies 
how frequently the product form basis is to be reinverted. Each 
time the basis is reinverted while the BCA portion of the code is 
in effect, Gl, G2, and BA, BB are recomputed using the new 
basis inverse. The same thing could be accomplished by DD calls 
to SUPERB with KEY = 1, and (PDI, IS) appropriately specified. 
SUBROUTINE ENDPNT (JS, PD1, IS, NET): implements the path-follow- 
ing algorithm described in Section III.4 [2}. The initial point 
is in a facet of the cell defined by (6) and the bilinear inequa- 
lity; we call this the initial facet. The algorithm follows 
the curve defined by the DD] bilinear constraints and the DD 
superbasic variables. The objective is to follow the curve until 
it intersects another facet of the cell and to find the point 
where that intersection occurs. We call that point the endpoint 
which is contained in the final facet. The initial facet is 
determined by the equation described by (PD1, JS). The final 
facet is determined by (PDl, IS). NET counts the number of 
Newton iterations which were required in the various subproblems 
of the algorithms referred to above. 

A common block is referred to for the first time in this sub- 


routine: COMMON/NEWT/H(10, 11), X(10) 2(10), ACC(2, 10). The 


following vectors are used in this routine: U(3), V1(10), 


(V2(10) F(10), DOT(4), RHS(10), UL(10, 10). 


H(10, 11) stores the jacobian for the various Newton subproblems 


X(10) contains the superbasic variables as X = (PI(MU), XX(NU)) 
Z(10) is the correction to X supplied by the Newton method, 
Z=H OF 


ACC(2, 10) sotres the two vectors which determine the current linear 


approximation to the curve. £1 (0) 


U(3) stores the coefficients for any quadratics that have to 
be solved 
v1i(10) holds the gradient of the functional defining the initial | 

facet 

v2(10) contains a vector in the null space of H, a tangent to | 

£ 1(0) i 

F(10) holds the values of f and any other functional involved | 


in the current subproblem 

DOT (4) a collection of accumulators 

RHS (10) a work space for finding the tangent vector by solving 
a linear system; RHS = H(-, ICX). 


UL(10, 10) stores the LU decomposition of the current jacobian matrix 


Subroutines called: QUADS, DSENT, CONCHK, FTN, DERIV NORM, 
DECOMP, SOLVE. 
It should be noted that the Newton's method being implemented 
is a modified Newton's method with descent. The iteration is 
k+1 x* 


ght a g® et tty | fap) ew) yy He 
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where a is the first scalar of {l, $s 4, > see} such that 


k+1 


fete yl Hecx®yt 


The basic loop is something like the following 


DO 100 K = i, 15 
CALL FIN(F, X) calculate f(X) 
CALL NORM (F, Sl, DD) set Sl: = I£(x)! 
IF (S1.LT.TOLCV) Go to 800 terminate if $l < 10> 
CALL DERIV calculate H: = f£'(X) 
CALL DECOMP(DD, H, UL) decompose H = UL 
CALL SOLVE(DD, UL, F, Z) solve by back substitution 
UL*Z = F 
ALPHA = 1.0 
CALL DSENT (ALPHA, S1) implements descent on the 
norm of f. 
pO 95 I = 1, DD 
X(1) = X(I) - ALPHA*Z(1) 
CONTINUE 
CONTINUE 
TAKE CORRECTIVE ACTION IF NEWTONS 
METHOD HAS NOT TERMINATED WITHIN 


15 ITERATIONS 


(10) 


(11) 


(12) 


SUBROUTINE QUADS (U, “MAG, ALPHA, BETA): finds both real roots 
of the quadratic 


ul) - a? + U(2) «a + (3) = 0 


If there are no real roots, the flag IMAG is set equal to 1. 
If there is a double root or if U(1) = 0, the root is stored 

in both ALPHA and BETA. If there are two real roots ALPHA 
stores the smaller root and BETA the larger. 

SUBROUTINE CONCHK (GMIN, KGMIN, MP2): checks whether the con- 
straints of (1) are satisfied at the point defined by the current 
superbasic values. The basic variables are evaluated in terms 

of the superbasic variables as suggested by (4). GMIN_ stores 
the minimum value of these variables and (MP2, KGMIN) indi- | 
cate which variable attains this minimum. The theory tells us 

that we only have to check the current bilinear inequality among 

the bilinear constraints described in (3). If this value is 

less than GMIN, MP2 is set equal to zero. 

SUBROUTINE FIN(F, X): evaluates the first DD1l bilinear 

functionals as they were defined in BLCONS (3). The last com 

ponent F(DD) stores a different function value depending upon 


which subproblem is being solved in ENDPNT. 


If MFLAG = 0, then F(DD) = 0. 


' 
th 
2 nl | 


(13) 


(14) 


(15) 


If MFLAG = 1, then F(DD) contains the value of the functional 
specified by (MPD, JJ). 

If MFLAG = 2, then F(DD) = X(JJ). 

If MFLAG = 3, then F(DD) = (ACC(2,°), X)> - BT (the normal 
hyperplane to ACC(2,-+), if F(DD) = 0). 

SUBROUTINE DERIV: computes the exact jacobian of the function 


calculated in f. The form of the jacobian is as follows 


PI (MU) XX (NU) PI (MU) XX (MU) 


D1 DIAG(PI (MU) ) *F1 DIAG(F1*XX(NU) + El) 0 


D2 -DIAG(F2*(PI(MU) + E2) DIAG(XX(NU) ) *F2 0 


The last row of H contains the gradient of the last functional 
of F as described above. 
SUBROUTINE DSENT (ALPHA, PNORM): Implements descent on the norm 
of f in Newton's method as described above. ALPHA is cut in 


K+l)q << pworm = We(x) or until ALPHA < 107. 


half until If(x 
SUBROUTINE DECOMP (NN, A, UL): (Borrowed from Forsythe [3]. 
The reader should look there for a complete description.) 
Implements a Gaussian elimination scheme with partial pivoting 


to produce a LU decomposition of the matrix A. A few statements 


were added to this subroutine so that KDET would change sign 


greene 


every time a row interchange was performed. KDET will eventually 


be the sign of det A (see Elken [1977a], sec. III.2). 
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(16) SUBROUTINE SOLVE (NN, UL, B, X): is also borrowed from Forsythe 


(3) with the additions of statements which change the sign of 


KDET every time UL(I, I) < 0. SOLVE uses backsubstitution 
to solve the system L*U*X = B. L and U are a lower and 
upper triangular matrix which are stored in the square matrix 
UL. 

At the end of SOLVE, KDET is equal to sgn(det A). e 
SUBROUTINE SING(IWHY): is also borrowed from [3]. It is called 
from DECOMP and SOLVE to print messages concerning the sin- 
gularity of A depending on the parameter MIWHY. 

SUBROUTINE DEBUG(MODE): prints information which may be useful 
in debugging the program when changes are made. The information 


which is printed depends upon the value of MODE. 


MODE = 1: Prints the time left in this job step by calling the 
system subroutine LEFTIA(TIME) where TIME is a 
REAL*4 variable. When running this program in WATFIV 
or on a system other than SLAC, this portion of code 
should be removed, 

The arrays BA, Gl, BB, and G2 will be printed by 
rows. 

The index sets defining the primal basis, JH and 
KINBAS will be printed. 

Information about the location and value of the first 
IH primal and dual variables will be given (KNU, KMU, 


NVH, MUR, €XX(I), PICT), I #1, «.2, TH)). 


MODE = 5: The coefficients of the bilinear constraints will be 
printed: BF1, BF2, Dl, D2, El, E2, Fl, and F2. 

MODE = 6: A variety of variables will be printed out which are 
of use in debugging ENDPNT: PD, KMU, KNU, INEQ, MUH, 


NUH, JH, DIGMA, PI and XxX. 


The remaining 1815 lines of BCA contain the subroutines which 
comprise LPMl. The reader is referred to [4], [5], and [6] for docu- 


mentation relating to them. 


II.4. Sample Problems 


The sample input and output are given below for two small pro- 
blems. The computational results for these problems are given in [2] 
Chapter V. 

For the reader to learn how to formulate an equilibrium problem 
so that it can be input and solved by BCA, (or HRA) it is necessary 
to read Sections II.2 and II.3 of [2]. We will discuss the definitions 


of D, b and C of (1) for a small problem here. 


The problem devised by Mas Colell (Wilson [7] is concerned 


with three consumers and two goods. By using the theory of Section 
I1.2 of Elken [2] we see that solving the equilibrium problem devised 


by Mas Colell is equivalent to solving (1) with 


-1 0 0 0 -0.9 
0 =] 0 0 -0.95 
d= | 0 0 -1 ot, b= | -3.92 : 
5 1. +25 te 3.0 
1. 5 -20 slg 3.0 
and 
-0.9 0 0 
0 -0.95 0 
c=] 0 0 -3.92 ° 
1 ue z 
1 de 1 


Notice that the sum of the columns of C is b. This will 
always be the case. We are almost ready to type the data cards for 
BCA, but we have to add a row to D to represent the objective row 
and assign names to the rows and columns. Below is our choice of 


names. 
ACT1 ACT2 ACT3 MEXP 


UTIL1 [-1 
UTIL2 -1 
UTIL3 “1 

~ G0oD1 5 a 25 Zs 
Goop2 1. 5 .20 1. 
OBJ t 
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We only enter the nonzero values to emphasize the fact that 
only the nonzero elements of D and b_ need be entered. 


In the sample input below (Figure 1) we specify that we want 


three consumers (TH 
two goods (L 


no production (IPROD 


C matrix input check (ICECHO 


Notice that we specify the Ith column of C by reading the Ith ele- 


ment from b and placing it in the corresponding position of C(°,I), 


and we read in the 4th and 5th elements from the card following 


Sample Output 


The output for the solution of the Mas-Colell problem by BCA 
is on the following pages. The current form of the program is rather 
verbose -- an option should be put in to print only the final solution 
if that is all that the user wants to see. 

The first page of output shows the values of all parameters, 
gives the problem name and the output concerning the solution of the 
auxilliary liuear program. Messages concerning time can be ignored 


in this printout because no timing routine was being called. 


OY Doi oO YHA owt - oS 
. 


—s Oe =e 
=“ 28 s «@ 


aq 


an 


MRS-COLELL 
&PRRM1 TH=3,Le2,TPRCC=0 .ICECKC=1, ENC 


6, SEnt 


ee aad 


SPARAM ICBG= 
HANE MRAS-COL 
ROWS 
E UTI 
L UTEL2 
L UTILS 
L 6cctt 
L Gccte 
N OBY 
COLUMNS 
RCT1 UTIL -1.0 acct o.S 
ACTI 6ccece 1.0 
RCT? UTILe Sta6 BCC 1.0 
RCT2 BCCCe2 0.5 
RCT3 UTILS sf. 0 BoCcct 0.25 
ACT 6oeccbe Q.20 
MEXP 6001 1.0 &CCCe2 1.2 
MEXP CBs 1.0 
RHS 
ENCCK UTILI =O UTIL2 =095 
ENDO UTILS “oe Fe 
ENCGH 60ct} 3.0 6CCb2 3.0 


ENDATR 
0 


a ee ee 


The following pages give information concerning the behavior 
of the bilinear complementarity algorithm, and finally, the solution. 
As a measure of the work done by the algorithm, the number of calls 
to subroutine ENDPNT, the number of scalar function evaluations, and 
the number of IX X IH jacobian evaluations are printed. The solution 


is printed in a common format for linear programs: 


JH(I) The names of the primal variables which may be non-zero 
in the equilibrium solution 

VALUE gives the value of these variables. The value of MEXP 
should always be near zero in an equilibrium. 

ROW NAMES is self-explanatory 

PECL} gives the value of the dual multiplier associated with 
that row. Those PI's associated with commodity balance 
rows are the usual "equilibrium prices." 


RHS bives the value of the original RHS vector (b). 
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Sampic Problem: _Whisman 


Next we give a sample input and output for a problem formulated 
by Al Whisman (Wilson [1976]). ‘The problem involves four traders and 
three goods. Each trader is concerned with two to four activities. 

To describe the problem we give only the D, b, and C matrices of 


problem (1): 


Again, only the non-zero elements of the matrices are given. 

The input follows. The first card, again, is just a header; 
it is not to be included in the input stream. 

The sample output for the Whisman problem also follows, page 


40. 
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WHISMAN 

SPARMT TREY Led, 19ROR=st ICECKG=t,. SENG 
SPARRM ICBJI=-S, SENG 

NAME WHISMAN 
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BeCC2 
UTILS 
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C3RCY sCccs 
CYRCI UTILY 
CYRCI BCCC2 
C4YRC2 UTLEY 
MEXP 500C1 
NEXP 80003 
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II.5. Source Listing for BCA 


The following pages are a listing of the program BCA as it 
is written in FORTRAN IV. The subroutines which make up LPM1l 


(Tomlin [1975]) are omitted. 
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$456 END 
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WE ARE GOING TN CALCULATE THE DOD OF COLFF 
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tOD= MESASCK i) 


Nos: 
an 


c 
FTreEar RAS ht 


OSuMr= y(Mu(T)) 


ANAS S-KC AC 
- 


9990; 
“ 


Cc) GO 1G 590 
UM+ Y(CJ)*G2(TDe!) 


Te ee ee ee et ee OL 
RROwCIntUlbuwnneocn 


FEELS EHH ES ELE 
core teo ere eee e eee 


NS 


70 


B0 


PI¢ 


AAAA 


120 


a 
Nfe 
wo 


~~ ee 
Wh 


Nr Our 


140 


T 


AANAD 


TF CICRALCS)0E929) GO TE 70 

RFI(K)= BFICK) + YC J)*FPECICBAS(I)) 
CONTINU 

IF (KLeNE*LD) GO TO 2&9 

INEG= 2 

Si b= GLi=t 

CONTINUS 


NOw THE StCCNO TYFE CF ECUATTION 15 CALCULATED BECAUSE 
K) IFS EASICe 


[FF (KNU.EG 2) PE TURN 
00D 145 K= 1eKNU 
K2= NUK) 


FE Che eCrelF} (CG TC 4s 

IND= [Ci AS(K2) 

HLe= BiLe+ 1 

Ir (KMU.EGe2) GO 10 122 

DO 120 T= 1LeKMU 

F20CK.1)= G2CI10Ne1) 

k2(K)= 1+ Acton) 

CALL UPAKC(K2) 

IF CKMUst Ce 8) GC TO 222 

9D 13272 i= te+kKMU 

D2 CKel)= Cet 

IF Q€MUE1) skEs My O2CKeT)I= YIMUCT)) 
10 125 J= Lem 

Te€ IDBAS( IV sEGEC) GO TO 125 
D2(KsL)= CACKeT) + YOU *GCACTICBAS(J).1) 
CONTINU! 

CONTINU! 

HFSC(KI= Oe 

DU 149 J= 1M 

IE CIDBAST IV ere 0) GG TO 249 
PFO(K)= AFECK)+ YC J) *BROCIDEAS(Y)) 
CINTINUE 

LG th 2s NE« OG) -CEe FC £4S 

INF G= 2 

Soe 

CUNTINGE 

RE TUN 

END 


SUBRCUTING &§ INDP (FP OLel So) 


AIS SUBPCLTINE CECOSES THF VARIABLE Yo SECOME IMPLICITLY 


BASIC AS THE (NG WET FEE FPEVST ECEKMERT GASGEST IN 
ARS 


OLUTE VALUE 
IMPLICIT] REAL *H (A-H,O-2) 
INTEGER FC PDL eP UP sSSeFR oe 7FLAG ORS eo Fe COPOO! -ELI eft 2 
INTEGER*2 JH( 3597.91 GMAC9S27) eK INBAS(1302?)¢ fORASCI 392) 
CUMMONZL NCCNS/GL ©4650 410) 656244600610) -FAC350),928(409) 
CUMMENZINDXIZ NUFC1LO) es MURCLO)eNUCIO) »sWUCLO) 
COMMORSIRNE XA SHel TGMAgKINGASs [ORAS 
CUMMENSUINS THN oe OD e KMUsKNUS BL 1 SL SO oMP 1 eNMe INOS IN 
COMMOENSTOLERZ THLE 2¢ TOLEDO s TOLEOVeSTEMXsSTPRD 
P= J 
EE GPO LVRCeHAt) Ge FG 36 
LOD= K IKLASE( TS) 
moO 20 I= LekNU 

52 


CABSSEGICICO e1)d 
(COMO eLECRIG) CO TU 20 
WEG= COME 
P= UCT) 
20 CUNTINU 
LE 4@O 36.4.7 ss THULFZ)? PF 
PE TURN 
Sy IPSs MC -ASE ES? 
ON 49 J= 1eKMU 
If {Pehi -G) GO FC 35 
sIG= DANVS(GPCIDD.J5)) 
P= MU J) 
35 COMF= HAPS (Ge(fOC.J)?) 
IF €COM@.LF«eBIG) GO TO 43 
BIG= Cims 
P= mMUCJ) 
40 CONTINU! 
{F (CBFC.LTeTOLFZ) P= 90 
RE TURN 
END 
SUBRCUTINE PIVOT (SSefF) 


THIS FOUTINE PERFCR¥S A PIVOT ON THE PRIMAL SUPER- 
BASIC COLUMNS [F POY.ZEQels OVAL IF PDeFQs-1e THE PIVaT 
RRINGS CCLUMN SS INTO THF PASTS AND COLUMN @F OUT OF THE 
BASIL Se 


IMPLICIT |} FAL*8 CA-hHsO-7) 

INTECER PO ePONLPLe SS eRR eZEFLAG ERS e Me ON PDD1 A VLLLOBLP 
INTEGEP* 2 JHO359D) «DI GHACSSE2) eK INFAS( 1307 )e IDF 4S(1392) 
INTEGER 2 LSTYPE eLAelL&« Ae IF ePUNSL CC 20) 6IC (S09) 
DOUALEF 2RECISTON EC secco) 

REAL A( 4907) eC CET) Fe CMI Ne COND cFEMAX ¢ SUMINE 


CUMMEN 1YSUM*eDOPROD OY e CE elP of (359 2 e XC ISO) oe VC 35E9)s 
LANE eCMIN s COND EE NAXs SUMINE « ICNAMC 1 302 02) eNAME (20) © 
INTE MP (27) eKINGO GL TIMeUTIMe I TINVeJTINVeMSTAT el O3Je [RIWP Oe IVING IT VOUT s 
STTONTe INVER OQ SLTTRLIM,» TFEF 2. JCOLPOeNROCW eNCOL eNELEMSNETASNLELEMNLEIETA,S 
&NGELEM eNINFeNUEL EMeNUFETAeNNE GOSeNL INES ISTYPE (359) 6 
DLAC1 292) sLEC2CO2 ) ePUNTE) « 

FIPUNC eNOE CT«eNDUAL NIP IT We IFEASs IFCRSOH 

CGOMMEN IToHe LICH AsTF PI WT» FFREG eK OUTA 
COMMON TA(40C0)- IECROOO) 
xx 
21) 
de 


ZB(359)- 
) 


CCMMON/LEL/PIC13C2) 302) 
CUMMON/ZLNCONS/G1I CVEO 


COMMOCN/SINOXIZ NUFCIO 


99 4G20400019) 66 350).98(490) 
10) 


AC 
MU ( 
aM 


COMMON INEX27 JH eDIGMAecKINEAS eS IOFA 

COAMONS SIMS THeN eM eb O eKMU ec KAUSPL Le OL? 

EPS= 16.%*(-13) 

RR= KINNAS( RR) 
t 


dF 
MURCIO) eNUTIO), 
G 


P1LeNMs INOs IND? 


TF CSS: 6(sGT es EH : 
IF (KRURCSS).NE6O 
CALL UNPACKISS) 
CALL FTRAR(1) 
25 €RBI= -—Ler 
90 TC I Lom 


ratte eee 


ANARD 


1 9 
20 


Ad 
82 


leu 


IF (CTeNELRE) ZECTI= YCTI*ZAC0RB) 


CONTINU: 

GO TG 3% 

JS= NuUF( SS) 

CACRB)= 1Le7G1 (RBIS) 

00 25 I= 2M 

TF Ob eN oh) LROETDI=H -ECLOTe IS) #Z3(RIED 
CONTINUE 

TF ({KNUCEC st) GO FO SS 
DUO SIO J= 1eKNU 

TF (NUCIJ) CEC eSS) GC TO £0 
V= GICF ed) 

GIORBe Js) = Ce 

30. 4¢C T= bem 

GLCTed)= €1¢01T2IJ)4 VeZE(t) 
GICSReJd= -GICRPP.J) 


CGONTIAU: 
Vz ACR) 
BACHRD= 9 


. 
OF €0 T= eM 
S4C1)= ACT) + V¥7ZROT) 
BACRi} = -EAC QR) 
WE An NOW PEVCTING CX THE NUAL SYSTEM, 
NOT AN IMFLICIT CASTC-TYPF PIVOT. WE CALL 
CALCULATE The PT VOT COLUMN. 


{FL= 9 

RH= [C0BPAS(ES) 

TF CRR «Gte EH} CE TE TS 
IF (MUFCRR).NELD) GO TO 99 
TF i= 

CALL SUMPEFECL e-1 oRhReO6C) 
JS= KMU 

Gu TO e@2 

JS= MUR( KE) 

ZA(RBI= Le /G20RA IS) 

6 6S t= Is 

IF (TeNt FE) ZAC 1)= -G2CTeISI*Z9(PR) 
CONTINU: 

IF (KMUeEC CD) -GO TO 116 
90 110 J= 1eKMU 

IF (MUCI) eEQeKR? GO TO 119 
V= GICRie J) 

G2ZECRReJV= Of 

DO 199 L= 1eN 

G2CTes)= Celted) + V¥ZE(1) 
GPICRBe J) = -GACPR YJ) 


CALL SUMERA 
PE TURN 
&ND 
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tf THIS 
SUTRA 


BRON nes PRREERE E soos 
Es —_— — 


a 


ANCHO 


aAnrexnaneen4 


SUBROUTIN® UPAKC (II) 

IMPLICIT EAL *8 (A-mHeO-2) 

INTECERS & JRO 350 De DIGMAC9IS2 eK INGASCL IOP VS INKBASTI IO) 
INTEGERS ESTYPF eLAsLt TAs IF PUN el C020) »1C(809) 
NOURLE MRECISTION E(85co) 

REAL A(4 


=C 
COO) CO S00) + CME N+ COND A RMAX e«SUMINE 
Ne 


CCMMOKR 8 SU DPROD OV e CEs DF oh (350 )e X( 350) 6VY (359) ¢YTEMP(350)- 
AcE eo CM iti e CONDO sERNWAXs SUMINE o LKCNAM( 1] 390202) eNAME (20) © 

aNTevec: PD) eSRINPSTTIMSITIMS ITINV eS TIAV eMSTATS IC Ise TROWP LIVING IT VOUTS 
SITCNTs INVER Qs LTRLIM, LFFE 2 oe JCOLPeNROW eNCUL wNELEMeNETA ONLI LE MeNLETA, 
4&NGELEM »NINF ¢eNUEL FMeNUFTAeNNE CD Js NLINESe ISTYME ( M59) 6 

SLAC L392) eL FC 2002 VePUN( ED» 

OSCE eee ene eee eee ee 

COMMON LICH, I ICHASTFPiWT se LEFREGeKOUTP 

CUMMEN TA ACACI0)s if (3000) 

COMMOCN/JHULCET/S&BF PCL) eRFLE LO) eE 14390) eEAC1 9) Colt C 


DO 100 I= I.NROW 
YCLI= Oe 
100 CONTINU 


Lie cecrEL} 
KK= ECCILCEY — 8 
DO 200 f= LieKK 
IR= ECtly 
Ue es Cre 

299 CONTINU: 
SETLRN 
ENO 


THE FOLLOWING SCUTINE MAKES THE CHANGTS IN THE INDEX 
SETS NECESSARY FVERY TIME A PASTS CHANGE IS MADF. 


SUBROLTINE BSCNC(S«P) 

INTEGERS 2 SHO 359 VD eDLCNACIS2AVeEKINFAS( 1302706 TORBASCLIIZ) 

CU4MUN/SINOXYZ JH eDIGMACKINR ASS TOPAS 

INTt GER © eSCAVARARSAV 

RSAV= KINFAS CR) 

JHOCRSAVI = S 

KINBRAS CH )= i@) 

KINGCAS(S)= RSAV 

SSAV= IDBPASCS) 

HIGWACESAV)I= R 

INBAS(S)= 9 

TOBASCH)= SSAV 

RETURN 

END 

THIS CURRCUTINE CAN CC THREE THINGS DETERMINED 8Y THE 
SARAME TER KEY", TT CAN ALT A CGLUMN,. REMOVE A TOLUMNe DR 
aOTH ALC ANC FEMOVE COLUMNS FCOM THE PREMAL OR DUAL SUPFR- 
FASTC CELUYNES. DEPFADENG CN WRETHER KEY (S le 2e OR Ge RESP= 
ECTIVELYVe f61 §3 1 IR =1 DEPENDING ON WHETHER A PRIMAL OF 
DUAL CULUMN TES SEING ADDED. PC2 [5S A SIMILAR FLAG FOR THE 
CUL UMM SBEING FEMLVEC e IS AND JS INDICATE THE PARTICULAR 
COLUMN TO F AOLED CR QEMCVED FROM THE GROUP HF CCLUMNS 
SPSCIFIED “1¥ POL ANC PRS. 


SUCROLIINE SUPFRE (KE VeFIIe Ef oPN20IS) 
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IMPLICIT VEAL A (A-H-O- 
INTEGER FI ePO1LePC2]¢SS 
INTECER* 2 JHO 350 eNLG 
Ack 

(2 

ds 


’ FLAG s8@SePel De 
INTFGER SC ISTYPE ol e 
: 
Cc 


4 
hed 

(S82 )e¢K INGA ASC 139° 
TAs TE oP UNL C1 20)6 
OOUPLE PRECISION € © 
REAL A(40C69).C0 F900 TNs COND e= PMAX »SUMI 
2P (359 ).X0350) 

e ICNAM(1 39202} 
TINVeJTENV @MSTAT 


COMMON SUM. IPFUD eDYs P 
ie 
I 
+e JCOL Po NRW pNCOL » NF 
N 
F 
t 


& 
A 
* 
Ce 
M 
F 
LAGE eCMIN es COND ee FL MAXe SUM 
2NTEMP (70) eKINPST TIM. ITIL 
JITCAT op LMVER Qs ITRLIM, TFF 
ANGOLEMs NINF eNUFLEM»NUES 
LACT SDP) eLEC ZCI? Ve PUNCE 
#TPUNRC »ROE GIT eNOUAL ».NIOTW 
COMMUN TECR,TTCHAsLFOIW 
COMMON LAC4S 
COAVENSL P12 
COAMON/SLENCON 
x1 
x2 
vf 


* 
! 
mM 
Eg 
A NEGOIsNLINESe¢ ISTY® 
) 

° 

r 


BAS.IFCRS 
FAEGeKOUTR 


) 
. 
gS 
A 
) 
iN 
0 
N 
. 
r4 
* 
. 
I 
. 


COMMONZINGD 
COMMONZING 

COMMEN/JOIN 

CPS= 16.94 = 

IF (KEY.FQel Ty 50 
IF CPD2.GT.) 

KM= KMU 
KMU= K+ 
KE C25 «4 
MH= MUH ( 
ee ig ee 


- Te OMA 
wee 


NUHCIM)- 1 
oN 
G2CJel+1) 


ar EREO t- S¥ 


RE TURN 
GO 
KN= 
KNU= - 
KF €3S «6 
NH= NUF(C J 
NUH(JS) C 
IF (NH ee! QeKN) GO TC 26 
00 25 [= RRekKNU 
JN= KRUCT4#1) 
NUCT)= JN 
NUH(CIJN) = RUHCUN) = 1 
00 25 J: Lem 
GIC Se T= ClOde IT +1) 
CUNTINU! 
NUCKN)= 0 
DO 30 I= 9M 
GLICTsKN)= 020 
IF (KEYVeECe2) 


1 

Te tH) CO TO 2A 
$3 

c 


FF TURN 


pid gb Aap 
MeNLETAs 


7086 


Ve OIGnrrrave 


~ 
camcd 


MO PG RG 0 me me ee ee ee ee 


NNSVNNN SSN SNS 
Uf 


N 
. 


Fa? 


1446 
745.6 
7406 
747. 
7436 
1446 
75 2d6 
ole 
tS2 se 
7536 
1546 
Fae 
Tae 
757.6 
(Re 
T59¢ 
(69.6 
“Ole 
Tue © 
TH36 
Toe 
1658 
106 
767 6 


[ok a) 


aan 


afan 


ovate latake! 
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is) 


70 


iss, 


ES 


an 


990 
132 


13 

-¢c! 

TH 
119 


139 


NOw WE ACC Yeo SUPFREASIC VARTASLE SPECIFIED BY POLeISe 


IF (PDO1etLT40) GO TC 7E 
KN= KNL 

KNU= KAU 1 

{© €ESsST«~tH) GO-TO ¢¢ 
NUCKAU}= F656 

NUHCIS)= KAU 

CALL UNVACKLIS) 

CALL FIV ANCY) 

DO of F> 1M 
GI(CIskKNU)= -¥O1) 

WE TUIN 


NOw wr AGRE BREIAGING & CUAL SUPERSASCIC COLUMN INe 


KMU= KMIi+ } 

IF CESsG cI GO TR F3 
MUCKMU)= 1S 

MUH(IS)= «KMU 

IF (IS eLF eM) GO TC 110 


THE ENTr&R ING VARTASLE IS A ZETA TRE ENTEP ING COLUMN 
A RCw CF INVOA) AKO A FCW OF ENV) #06 


IDD= KINEASCIS) 

O01 75 $= 1eM 

YCL)= Ge 

YC TECY= Le 

CALL BTIVAN 

CALL SRIF TRU 3227) 

90 &) T= 1. 

te CLIDE4&SCT2.EQs.9) GE TA ad 
G2C INRA (C1) eKNU)D= XCOT) 


CONTINU. 

Offs 100 J= MP1 eNM 

fe (CPOCAS CI) «EQ. OG} GO TG LEC 
CALL URMACKEO I) 

NOT= Jed 

3 SO I= 14M 

If CICEASGCIF2+C960) GO TC 90 
OUT= DCT# VCEIIeXECL) 
CONTINUES 

G2CILSASCIVeKMUD= OCT 
CUNTING 

SE FURN 


THE ENTERING VAP TABLE [5S A Pte THE ENTERING COLUMN 
PARTLY CI#INV(IE) AND PARTLY C2+ CI*INVIR) €D. 
“INV(Q'3) TS PART CE THE GASES INVARSE COFRE SPONDING TO 
K-3AS TC COLUMNS AND THE TPASIC ROWS. 


ToO= KIN 4SC15) 

OO 149 T= Lem 

YCIT)= Ge 

Y(rCod= 1. 

CALL BIPAN 

CALL S®TFTR( 22?) 

nM 14) 1 1.” 

TF CICEAS CE) e#V09) SGC WO 146° 


5? 


7AG6 G2CLTORAS CT) sKMU) = XT) 

TOS 149 CONTINU! 

7706 90 160 J= MP1.NM 

771 TF (ICEAS(C I) eEQe%) GEC TO 1&0 

7726 CALL URVACK CS) 

7736 OOT= 0.% 
$14.6 bO 199 t= Fem 
q lt eo If (ICKRASCI)-EQ-C? GO TC 180 
76. DOT= DCT+# YC T)*xX (1) 
: Pt ay 1590 CONTIAUS 
7786 G2C1OBAS( J) eKMU) = YC IS)+ OCT 

779.6 1660 CONTING! 

7AO” RETURRA 

7TBle FNOD 

7H. SUBROUTINE RECALC 

733. C 

TH4as € THIS SUFROUTINE RECALCULATES TRE SUPFERRASIC CCLUMNS 

1596 ¢ USING TRE CURRENT BASIS IN CR FORM. IT CAN 4L5C AND 

THEs € A VARTABLE AND COLUMN TC THE SUPEPBASICS. 

7E7 « Cc ‘ 
| Ft. IMPLICIT REAL*B (A-H»C-2Z) 
7B%° INTE CEE PN SPD 1 sO C2 .sSeRE eZFLAG eRe MeL IDM BLL AL? 
79C © INTEGER? JHC 350 De CI CMACSS 2D eKINBASC 12027) 6 1DRAS( 1302) = 
) 791. INTEGEP X22 ISTYPE sLAsLlE elAe LF ePUNGL COC 29) 26101409) : 
S426 DOUELE PRECISICGON ££ (Ff 9C0) 
79S PE AL £04099) 6C€C £99) eo CMENSCOND cE RMAX o SUMINE } 
; 75946 € | 
7956 COMMON DSUMsDPROD DY eDE OP e2C3S I) oe XC FEDIV CVI ISI VY TEMP (359), | 

Vibe LAs& oO MINe CONC eERNAXs SUMINE 6 LOCNAMC1L 30 Se 20 eNAME( 20) 

7976 ANTEMP C20) KINPeL TIMeITIMeLTINVeJTIANV sMSIAT SIONS es FROWPCIVING IVT S } 

(IR. STTCNT e INVERGs [TRL My LFFEZ » JCOLP WN ICW oNCHL @NELEMeNETAe NLELEMeNLETAS 

199% © 4NGELEM «eNINFeNUEL EM eNUETA eNNE GIO NL INE SS ISTYPF C359) «© 

aIGe GLACT3O2) LEC 2002 ) PUNTA) « 

a0le 6 ITPUNC eRUECITe+eNOUAL eNIP TWe IFEASe IFCRSH 

BOD. CCANON TTCHe LT ICHA,IFF I WI e fSFREG eKOUTA 

BI3- COMMEN TAC4099)6¢ 16 (4790) 

n246 COMMOERSLPLSPI CL 3C2) oe XXK(1 792) 

aISe COAMONZENCONS/ZSL CAS 3010) eGPC4ID01 I) EAL 25%) oF 3(490) 

lerery COMMOR/INDXYZ NUECIC) eMURCIC) + NUC(L9) »MUCLO)D 

OW. COMMONZINDX2Z JH eDIGMAeKINAAS*s TIDBAS 

rh Jtte CUAMINS LOCALS BT oe Ne SIS eo MELAG DD 1 9 M0 1) eMPDG INE Qe KFIING K JAC 

4996 COMMEN/SOIMZ THeNeMeDD oe KMU se KNUSAL IT oil 2eMPisNMe INOS IN 

416 EPS= 16.**(=-13) 

Hille It (KNUseEC eI) GO TO 26 

Mlee OG 28 v= LeKKU 

Hi je LodD= NUT) 

mola. CALL UNPACKITECD) 

nlSe CALL FT!?AN(1) 

ole. OO 25 [= Lem 

t17 25 Gitlesi=2 = FC? 

alBe 28 CONTINGSE 

Ai%.% 2c CALL SRIF TRO1¢63) 

b20 6 CALL FTWAN(1) 

APL. pO 22 f= 1M 

Oo Pe Je CRE TPS YC 

Aa’Se 3) LF (KMULEGeD) GO TO 149 

te NY 2279 J= Lek 


O25 InpD= MUCJ) 
AM fF CICHeLE ow) GO TO 7S 


nan 


40 


60 


70 


30 


190 


16) 


ZETA= VARIABLE IS 


1aM 


Fe) 


2&929) 


{F CLEGAS 
( Jd= x 


G2CICBAS 
CONTINU: 
90 G2 K= MPI 
fF CECRAS (CK) 
CALL UNPACK( 
NGT= Oe 

bo SO f= le 
iF CIChacc! 
HOT= OCT+ ¥ 
CONTIAU 


NM 
FQ.) 
) 


M 
».U%. 
CL) *® 


G2(1OBASCK)I.eJ)I= LO 


CONT TNUF 
GO TO 1 30 


0) 
x 
10] 


ENTERING 


co TC 4 
1) 


GO ¥O 6 


GG TC § 
) 


PI= VARIAELF [5S ENTERING 


Q F= Lem 


AS(K) 26900) 


ro) 

§ 

ee 

a, 

6 

At 

EO Ka MPL aNM 
Of 

L URPACK(K) 

y Oe 


fF CIREASCTI) eEQ20) 
ONT= vwCT+ Y¥OT) exe 
CUNT INU 
G2CLOR AT CK) JIE YE 
CUNTINGE 

CCNTINU~ 

DU 159 T= 1M 
¥EL}2 O60 


Tr CK IN ASOKM) NES 
CALL AT’AKR 

CALL SRIF TPC 32?) 
it} BOC T= Let 


IF (CIDUASCI) «NE eC) 
OD 179 J= WPT eRM 
IF CIDPASTI) cE Qe 2) 
CALL UNRCACKOD) 


DUTS Oo 
DO 363 f= Is 
TF (1OBAS ECL) cE C) 


GO TC 120 


9 


9 


iol 
4 


1oC)4 oct 


Cy VYOCK TNPASONM) J = 


BSCIOSAS(T))= 
GO tO 176 


GE TC 165 


aSAe NOT= DCT+ XCITI)*v C1) 
B39. 163 CONTINU: 
A906 BACICBAS(I))b= DOT 
AIL. 170 COUNTING: 
BOP. RETURN 
A936 END 
mG4. SUBROUTTIAE ENOPNTOCISeOPC1l eI SeNFT) 
BID. IMPLICIT? REAL*8 (A-HeO0-Z) 
BI. REAL*3 NINeMIN2 
mile INTEGER +2 I[STYPE sLAsLE stAs LE sPUNeL © ( 2036 IC(H0)) 
a93. INTEGER*&®2 JH( 359) sDIGNA(GS2) eK INSA501307),I9R AGC 1 392 
OID. INTEGTR PO ePDIL PCS oe SE ekRR oe ZFLAG ei. COLN1 BL Ie Fite 
330. REAL C300) 
YOle COMMON/NEWTS HOLC e112) 9X610) 67010) .ACC( 3010) 
Wee COMMOEN/LEIZE TCL 3IC2Z Po XKO1L292) 
VO3. CUOMMONSSL CST SBF 1 (LO) eEFT CLO) eCELICLOD cE2CLO) eC of Cot C 
W9O4-e COMMONSELCST2ZDI CIC oe lLOeEA(GeINAD oF LC Sel MdeF 209019) 
40S. COMMEN/ZLNCCOCNSZGL O75 06199 eGF(49I019) FAL 359) 6 44 (49)) 
9066 COMMONS IROKIZ NUEC1O) MUFC IODSNUCIO) »MUCTY) 
VOT. COMMON/ IND X24 JHsOLGMA.LKINEAS, LOBAS 
9OKe CUMMONSSCALS ATe RR eID eo MELAG CDN Le Pe PU eMPDg INFOQ yg KF UNS KSAT 
WOOD COMMONS DIM THeNe Meh De KMU KAU SSL Ie OLA eM 1L eNMe INQ. INP 
Gide COMMPON/INTZ 19S C27 +e KCET eKOCUNT eI SING 
Y1lle COAMERSTOLERZ TOL Fe TCL CC e+e TOCLCVeSTEMX,STORD 
G12. DIMENSION UCSIeVLCLIMVECIC) oF C10) sDOT( 4S) e2HST1C) eUL (10010) 
Yilse € 

; G1i4e € THIS PROGRAM WILL FIND TRE ENDPOINT NF A CURVE 

| ise Cc DEFINED BY A SYSTEM OF DC-1 RILINFAR FQCUATIONS IN DO 
GVlGe Cc UNKNGUWNS e CNE ENDPCINT FS KNCWNs AND TRE OTHER END= 
917. ia PCINT [SF DOXLTFPMINED 3Y TEE FIRST INTERSECTION NF THE 
VIiARe € CURVE wWITek Chet CF NM41 ENEQUALITY CONSTRAINTS. 
41U. Cc THE SYSTEM ES 
G20 6 € 
tele Cc BE b+ OCEXCME)-— DIAGOX EMU) *® OF TX XOENUDF E1)= 0 
“226 C BFE+ DP*X CMU - DIAGECX CMY) D* CE 2EXKOMY D+ E2j= ie) 
3236 Cc 
W246 Cc WHERE CNE CF TRE AMOVE IS AN INEQUALITY IF ZFLAGeNF el Z 
G25e6 ¢€ THE INEGUALITY 1S TE LAST IN THE GROUP CEFINED AY INEQ. 
Wee Cc Tec LEINEAR INFQUALT TIES Apt DEFINED ‘11 LCWwe 
1276 € 
Y2B6e NETS 0 
YPDe FRACHLe 
Gade Cc MFLAG=0 WREN WE ARE S0GL VIG FOR THE INITIAL 3 POINTS 
G3le ¢C MFLAG=1 WHEN CRE SF THE G-CONSTRAIANTS IS RINCING. 
WIL 6 Cc MFLAG=2 WHEN OAS OF THE VARTARLES GOES TO ZERO. 
Q34e Cc MFLAG=3 weKEN WE USE TKE ACRMAL HYPESPLANS TO AFFINE THE 
G34 « € SUBPROELEM SEURT CF THE '80UNDARY OF THE CELLe 
V3Se ITER= © 
W356 1sGcT= © 
W376 MEL AG=C0 
9 IR oe L=1 
Y3%e KK= 0 
W406. 9bi= DD- 1 
G41 1O= ING 
W460 KMUL= KV¥UF IL 
G436 IF CINEG e&Qe 2) o= INP KMU 
G44. IF (FD ef GeC) ID= OD | 
4456 1¢x= If | 
4b t© CKMULFG 29} GO TH 615 | 
“W47Te DO 610 I= LeKWU 


60 ‘ 


; By ne anal oad 


1) 


34 


XOT)= FLOMUCT))D 
IF (CKNU.EC.O) GO TQ 625 
DO G29 [= IsKNU 
KM= KML¢+ I 
KCKM)> XXCKRUCT)?D 
ALPHA= -t. 
LTER= FTER+ 2 
IF (Oe s:T21) GA 10 6 
KNET= ft 
GO TeeAT 
TF (KNULEE 2D eARs KML.FGeDD) GO TO £30 
CALL NEVIYV 
G 


ISGCT= tSEEt4 1 
TF CISGCT.GFeCD) GO TC 359¢ 
OT 36 Fe" {sOO¥ 

RhSCLY= -HCT eI CX) 
CUNTINU: 


fF €Te€x.EC.00) GE Te 25 
oO -4 J TCx.01 
4Fi= J+ 4 
eCOl 
= H(1eJP1) 
CONTINU! 
IF €9991 .6T6 
IF (vuABS (HO 


i} GO Fo as 
ls 

ISING= 1 
Ff 


ib) «GE TOUCV) GO Yo 26 


VECLI= PHSO1)7HO101) 
fe UHeLGMISCEE. Ok COLTe 3S 


mOoiT= 2 
GO VO 4D 
CALL OLCTOVMECOOLs te UL) 
IF tiSinGsAE et} €O TO 2A 
IF €FCX wEOe 13 GE FO SE 
icy= CxX- 1 
GQ TC 7 
cx be 
GO fG 7 
CALL STIVEQ(PO1L.UL+¢eRHSeV?) 
9a 39 1 1,001 
ivé= C=. ¥ 
TF C8 TeICXI GO TC 41 
TF1I= IM+ 1 
vetlPld= VectM) 
CONTINUE 
Gu TO 41 
KDE T= 1 


DO €24 {f= I2DD 
V2ETI=H On 


V2CICX}= 16% 
lf CETF+eETet) GC TO 48 


IF €POe'G.0) GU TC 120 
te (PUslGe—-1} GO TO 119 
JE= KINDASCIS) 
IF (J%eNF oD) GO TO 194 
TF (VECTDO6GIeC) GO TN 15 


” 


GG TH 140 
DO 105 I= leKAU 
Ik= KMUt 1 
SUMV= SUM# V2C IK) #E1 CIM eI) 
CONTIAU 
IF (SUM.LT2929) CF TE 140 
GO tG 13= 
JB= FOREASOIS) 
fF {JReNEe«D) GO. FC 842 
if (V2¢T0)-.GTeI) GO TO 135 
Gt Ta 1490 
09 115 [= ieKMU 
SU%> SUM+# VeECI) *C2(JR6T) 
CONTINU® 
TF (SUM tT. C60) GC FE 140 
1 13S 


CAL CULAT®~ TRE GRADIENT CF TRE SURFACE se FC(OD = Dele 


TF CINE GeFQed 
FELO}= CE NGo ed FOTCMUD- OITCING)® 
INGeoe )*XXCNL)+ FICING) De 
INO= MUHCND) 
00 122 1+KMU 
C1CIANGel) 


SUM# FICINGC +s J) *¥X(JUK 
CING)- SUM- £EICTNO 


) 
) 


LuU3G.e f 

1049. 2 PIN®FICIKGQs, J) 

1041. 

134?. 

104+. = Sum SUM+ VICI) ¥*¥V2CT) 

1044.6 (SUM eLT2)) GO TO 149 

1104S. 7 Fe 235 

1046. c 

1347 « FADD) = BFACING) € D20KNUs e DFP TOMI) XC TPO (F2CINP ee )HPY (MU) +E 2) 
1048. 

1049. I8= INP+ KMU 

10506 90 127 J= IeKMU 

l1u51. SumM= SUM+ FECINPsJV*xX0I) 

1052. PID= XC1F) 

10536 DO 129% J= 1eKMU 

1)54. VICIJI=R DECTROPSII— FID *F2CINEs J) 

1055-6 CONTINU 

105. DU 129 J= KMU1eND 

1057. 29 VICJII= O20 

109%. VIC IS)= -SUM-E2¢C INP) 

10596 

1050. le 
lord é UN" su 
1062 f 20 
loos 

1064 p 

106% © INIT= -#DET 
lub6é. > CALL NERM( V2 
1067. Si= (KCETS* IN 


oo 
M4 V1ICT)I¥*¥v2¢T) 
»’ GO TC 140 


gape rene cer gn 


lubse 
L1u69. 
1LO706¢ 
LO7Le 47 CONTIN 
LW7]e E 5 
1U7 46 : ; CCACCC Ls I) 6 1 =203) eJH1Le DD) 
1074-6 = 
1075. : FINC TRE FIRST POUNCARY THAT THE TANGENT. OICALPHA) 
LU7bH. HITS. we APF TRYING TC FIND 
1077. THE SMALL SSF AL? tA SUCH TEAT 
LAA " GOL ee) FOQCALFHAD 4 LECT) = 
lG?7Xe 
10%Ce IFLAG= O 
1081-6 IF (CKNUSCEG eS) GO TC 249 
LI42 6 
lu3s3e 
LOS4G.s ef a IVeLEsiv) CO TO 250 
L065. 2G: 2% 
10AG6. 
1047. 
LUMH. 
10A%.° UCttdt+ GICL eK )*ACC(LL eK) 
1090-6 
1091. CORT ENUF 
10926 UCIy< ) 
109.36 ) FeatCerFZz? GO TO 250 
1O74,. 3 ) 
10% 56 (ALPHA : =TaUeZ3 Go ™T 259 
LO sre (ALPHA CLaO) GC TO 244 
LGU. 
109%. 
1499, 
TOD= TCO+ G11 eK) *ACC( 7. KM) 
TF (TOD eGTe -TOLFZ) GO TO 250 
If €LECAGLEGe 15 GG FC 243 
JFLAG = } 
JJ= i! 
MIN= ALPHA 
VPD- 1 
6C.t0 269 
fF CALPHA.GE eMIN) GC TC 250 
JJ= 1 
MIN ALPHA 
MPC= 1 
CONTIAL: 
IF (KMUceEC eI) GU IG 
Of) 259 T= feN 
IF (DICMACT) = GG TO 259 
BO 286 Lis ei63 


IN DD 
. 


UCLLI+ GACT sKI*FACCILEL KK) 


CONT ENU' 

UC2)= 

IF (LABSCU eTOLFZ) GO TO 2569 
ALPitA= 

{fF CALPHA CLFZ) GO TO 259 
TE CALPHA, ecy ct TO 258 


WHI DNF ewveo' 


crs 


9 
J 
) 
| 
J 
Q 
) 
1 
1 
1 
i 
1 
1 
i 
1 
1 
1 
2 
Pd 
4 
2 
? 
? 
<M : 


= 


| 1128.6 be PSPS Kis PeKRG 
| 11296 TOD= TOOt ACCLAcKI*G?CT.K) 
| 1130-6 252 CONTINUE 
| Lisle LF €TFTON «GTe -TOLFZ GO TO 259 
8132 6 298 IF (CLFLAGeFQe1)} GO 10 263 
t ELIASS IFLAG = 1 
1134. JJ i 
Li 7Se vit AL MRA 
Lisbe MP -1 
LISe GO TE 259 
1138.6 253 If (CALFHA.GE eMIN) Ge TC 2689 
LI3Vs JJ I 
1140.6 MIA ALPHA 
Liddle we =1 
L142. »9 CUNTISA 
1144.6 “ “MeL AC 
1144, IMAG 
Li45e6 red Lea 
Ll4o00 558 i) (f1)= 0420 
Ll4a7. ’ IF (ff JeECe2) GC TO Ott 
11456 0 ] f= 1,90 
11496 IF (14eGTs*KMU) GO TO S91 
115%. IOTOCLI= OOTCAP+ CLLINGsLIY*FACC(S.f) 
Rh -« OT (2) NOT(2)+#+ OLCINC sf) FACC(P6T) 
11526 Gt t¢ 892 
Libstae S35 Ir I- KMU 
11546 OCT(2)= DNTC S24 FLOINGs TKISACCT( 241) 
1155.6 DGT( 4) OUT C4Y+ FICINCeIKI*ACC( 2.1) 
115606 S32 CONTINU 
bist UCL )= -ACC(2,.INQ)*NDT( 4) 
11566 UCFI= =~CACECAsSING) FEL CINGI4 ACCOSe INCI¥IOT (4S PDF DOTI2) 
LELEVe U(3)= BFICING@I+ CLOTCLI= ACCC 3s INO) e(CCT( 3) 4 FLCENO)) 
LibQe CALL GUADS(U¢ IMAG. ALFHA st TA) 
L11Ale IF CAL PHA cL Te TCL CV} ALPHA= EFFTA 
LlH2e GO 7G EO 
LlOd te 6090 16> INP+ KMU 
LlOGe IF (KMY.EC.S} GO TH. EOS 
L1lOS.e DO €32 1= 1eKMU 
11566 JCVCL} COT Ctd¢€ CACINF of P¥ACCE3Sef) 
1iG7e IBTL|Ss= ONT (2)+ CECENP.E) *ACC(2 61) 
LimSe DET(S)= ODOT CI) F2t INP eT FACT %61) 
Li te DOT(4) NOT (4)4 (INF ae LD *FACT( 2410) 
117906 6902 CONTINU! 
Liv te 692 UCLI= ~ACCC2]e¢ ffl) eHOTI4) 
LivZe ud?) -(ACC(2esIA)*E2C INP) ACCC 4A.T NP) RDITC4) + DGTP?) 
be Roger UC3Z)= PF2ATINP)€ CETCLIN=- ACC (Se INGE €(COTC3)+ F2ACTINP)) 
1174. CALL GLACE (UZINMAC, ALPRALRFE TA) 
LET Ss IF (ALPHA eCLTeTOLCV) ALFHA TA 
Li 76.6 > IF CIMAGeEGel) G TO 26 
li! 77e LF (SETA e«L Te ~TCLECV) GO TG 260 
1178. IF (ALFHA el Te —TOLCV) ALPHA= RETA 
Li7D« [F (ALFHA eGTe MIN) GC TC 2E9 
11 de SIZE Mt | 
Lthl«e MIN ALPHA 
1182e MPED= } 
LTidte ¢ 
Lisa. ‘ NOw CHECK ¥0 SEE LF CRE EF TRE VANRTAELES SGors TO Z2Feo 
1185.6 ¢ BEFORE ANY CF TRE CCNSTRAINTS BECOME TIGHT. 
Liste é 
LIS7.s 260 MFLAG= 3 


DO 280 

pn 2ES ace 

(tL I= ACCKLL ef) 

CONTING 

IF €DABS(&G £23) «hTe FALCFZ) GO TO 2959 

ALPHA= -U(3)7U(? ) 

If (ALPHA UT. MIN eORe ALPHA et Te TOLFZ) GU TO 28) 


eo rvc cra 
NPunreCoc 


~ oF 


I 
“ 


= 
~f°9 


LP 


soIN 


15) MENe MPCs. JJ 
ETEPLENGTH= "ef 150e65" CORSTRAINY TYPES( tel Seta *soTSe*)*} 
sh. te StPORi “Gl TG RES 
MEK&® STPERD 
FLAG= 3 
WRITE (6eE3E) MINe MPD. JJ 
BT= 926% 
OU 290 1= 1-.bD0 
X€J)= MIR * ACCC] 
GT= GF + XI *XACC 
CONTINUE 


NN ee eR eee ee ee ee 
5 


Ovtls 


eT) + ACC(361) 
(2.1) 


NEWTON'S METROD FO? TRE PRORILEM OF FINDING THE INTFRSFCTICN 
THE CURVE WITH THRE COASTRAINT OEFINED SY MPD AND JJe 


KK= KK + } 
TF (KK .GFE e5) STNF 
OO 100 k= Ie1¢ 
KML= K- 1 
IkwOIM= OO1 
LF CMFLAG eG TRwOIM= DO 
CALL 
CALL NORMCE e 
wel T (6 45¢ vi D 
FOOMAT 2A ¢,9F19.5) 
weit (60936) Sl» worM) 
$F CSP Ae LOCC) 
CALL CERIV 
NFi= NES 4 
LF (SECT et 5 s$ 
IF (DABS GO TO §& 
EELS Fat 
( Te é¢ 
ECLMP CT RWDOIM oto UL) 
SINGeF GeO) GC TO 69 
(6594C) 


If 


STi: 
CAI SCLVEC IRWDOITM eULsF eZ) 
ALPKH= 1.29 
CALL CSENTCALPH.S1) 
IF (ALPRLGT. TOLFZ7Z) 3C TO 70 
WRETE €6eS94°) 
Ge T2 BOO 
ce <5 {i= }st0 
REL P= KROUVI=— SP 8 7ETY 
CCNT INUVE 
CONTINU 
MFLAG= 7% 
MIN= MINS e 
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SER Ky 


ANDAIAAG 


TF (MIRAGTeTOLHD) CC TC 255 
Gt TC 3¢C 
800 IF (MFLAG.EFQ.3) CC TO § 
CALL CUNCEKCGMIN sKGMIN MP2) 
[F (KMU.6C20) GC TO 146 
90 145 JT= 1:«MU 


KM= MLCT) * 
PIC(KM)= xXC1) 
TE (X01) -GE eC MIND GC TO 145 
GMIN= 2ACT) 
KGMIN=S KM 
MPe= -] - 
14 CONTINU! 
14€ IF (KNL-ECeI) GO IC 149 
Oo 147 I= KMLIEDE 


Ik= I- «MU 
KM= ALTIK) 7 
XxMEKM)= XC) 
IF €XC1) eGEeGMIA) GF TO 147 
CMIN= XE) 
KGMIK= KM 
ve?= ft T 
147 CUNTINUE 


14@ IF (@GMIN et Te -—TCLBO) GG Te 170 
EF (#PE okks OF GH FO ESS 
POL= 0 
RETURN 

153 £F (MFLAG.ED«2) CO TR 180 
IF (MPDeEGe -2) CO TN 1251 
1S=-JHCS JS} 

PHL= 1 
RETURN 

PSI FS= HICMALII2 
PHi= -1 
RETURN 


1S0 IF (JJect.KMU) GC TO 1£2 
S= MUGJI) 
PDi= -1 
PETUPN 
{S2 JJ= JI=— *MY 
1S= NUCII) 
POl= 1 
RETURN 


IF THF CURRENT VALUE CF X VIQLATES THE CONSTRAINT OFF INES 
AY KGNYIN AND MP2, FINC A GCCOND STARTING FOINT FOR NEWTON'S 
MC THOG &Y FYINCEIKG TREE FCIAT CN TRE LINE SFGMENT 
ACCEL ee) + ALFHA*(X — ACC loeel)de O<= ALPHAC= 1 
THAT SATISFE TES THE VIGLATER CONSTRAINT EXACTLY. 


173 JJ= KGWIEIN 
MPD= MED 
“FLAGS 1 
with Té (4 017%) MPC eS e GMINA 
174 FUPMAT(Z.* THE QLADHATIC PICKED THE WRONG CONSTEAINT se § ef 
3° TRE WCST ENFEASIBLE CONRSTRAINE EFS FYPF * 51 Fs 
»® 


NOMBE RS ef bes4o* WITh A VALUF OF ® fF 146m) 
Oy U7) f= eB 
i71 VGGCCV= KROLCE= ACCOSe £4 
0D i7 { 144 


66 


= 


= 


. 


ed ee ee! 


ee ee ee ee pe te ee et ee ee eee et ee et pe 


$14.6 


ADAH 


ees 


176 


DCICII= 9029 
IF (MPC.EC.O) GR TC 190 
IF (MPCsEGe-1) GE Yo bee 
DD 172 T= LeKNU 

Ik I+ Kau 

BETO E = por: £r1¢€J3J 

BEaTL2}= e 
CONTINU: 
STR= (-3406455)-00T(1)97/00TC2) 
Oly LFS i ehd 
XCTI= ACC(ZeI}+ STR*VICL) 


GO TG 2 


NO 14a3 T= PeKMU 
DGTCLI= COYCL)+ GECISeI)*¥ACCO36T) 
OC Pls d= COT (C2)+ GPCIIsITP¥*¥VICT) 
CANTINGe 


STR= (=-FECIIJI-“DOTO1II7ACGIC2) 
bY 196 T= 1-D0C 

KXOCLIV= ACCOSel)4 SIR*VITI) 
GO TO 25 


THE BYLINEAR CONSTRAINT IS NOT SATESFIED AT Xe FIND 4 GOOD 


ETARTIANG PGINT BY SOLVIRE TRE CUADPATIC DFFINEO AY THE INTER=- 
SECTIGN GF THE SEGMENT OEFINED AROVE AND THE CONSTRAINT. 


193 


209 


20? 
P04 


<os 


30 


TF CINEOeEGeZ GC TO 42C0d 
BO 292 f= saSE 
IF ..€BsGTskML) GO TC LOY 
NGF C1)= CUOT(1)+ CICINGsT)#ACCC 3eT) 
DOFC2= ONT C2)+ CICINGel)*V1ICT) 
C Th 192 


tK= [- KMU 
DCT(C3)= HaUT(syt lide atta at dicen ES, 
FICIN 


DCT{(4)= DOT (4)+ Pe TKIEVICT)D 
CONTIAU 
UCLI= -VICINO)*€OCTC4) 
UC2Y= -CVLOING) €ELCINGI4+ ACCC3e¢ING) *COT(4))4¢ DOOTC2) 
UC3)= HEICING) + OOTOC1)— ACCI3.INO) *CCOTC3)#+ ELCING)) 
CALL GUACE (ths TMAGeSTHSEE TA) 
Ie (STK.-LT-TOLCV) STR= PETA 
DG 193 f= 2.66 
XCPEV= ACCO3s1)+ STR*VICT) 


GO TO ave 
IS= INP+ KMU 
TF € KMUCFG.0) GC TO 204 


OU 202 I= 1eKMU 
UU 8 DOT (1)+ C2CINP eI) *¥ACC(3eT) 
DET (2)= DOTIP2)+ CACINPe ED ¥VICT) 
DCT(2)= DOT (3)+ F2C INP.) FACC CST) 
DOT(4)= DOTC4)+ FACINPST) *®VICT) 


CUNTINUE 


UCT)= ~VILCTIE) *DO TC 4) 

UZ d= -EVICETA IEE SCINE D4 £0003. INP ENETCAY De NOTE2) 
UC3)= BFZCINP)* DETOCII=— ACCCSsINPI*(COT(3)¢ F2CINP)) 
CALL QUACE (Us IMAG.STReGETA) 

LF CSTP.LTs<TOLEV) STR= BETA 


90 203 t= 1.00 

XCEY= ACCS e+ STR*VICT) 
G9 Te 2o5 
WRITE €559280) 


Labs &SO FURMATCIXe* THE CLRRENT LINEAR APEKCKIMATION TQ THE CURVY ISce¢ 


L469 Ve LefolXe* A2(*U)+ Aa") 

LlilJde 8S5S FORMAT (1X2! 17-7) 

C27 i< 94F FURMATCLXs*NORMCE CKD P= Sef 1Oefe® F= *,9F10.6) 

lsal2e Q4O FORMAT (UX eSS9** FKALGCRITHM PONHED WITH SINGULAP JACOKBITANS®* ee #) 
Sits. G45 FORMAT (1K e*NG RESCENT FOSSTELEF --— ALPHA WAS 7FEROD.¢) 
Late 9H0 FURMATCIXs "NEWTON'S METHOD FAILED TO CONVFR5F et) 
1475 e STae 

l.a7Ge ENO 

1477.6 ¢ QUALS fF INCS THE SMALLEST NONNEGATIVE POUT OF 

Li7be € U1 RAL OriARK? @€ LC 2) BALPFA ¢ UC3Z) = Oe IF THERE 
14796 € IS GNE« IF NOTs EMAG YS SET = Arc 

14306 € 

L691 « SUBROLTINE QUADS (UsIMAGs AL FHA,RETA) 

14327. IMPLICIT FFAL eR CA-t.G-2) 

1st45e INTFUCE® PPD DSPDP eS ePk e7FLAGePeNNe CO] eALIL EL? 
1.47346 INTEGE&Se® 2 F40.350 ) eDLEGNA(G52) eK [NFAS(C 130? )¢ FOR ASC 1397) 
L4H5e COMMCOCN/STOLEKS TOLFZ¢ TOLED s TCLOVeSTPEMXeSTPRYD 

144566 COCMBMCN/JD INS THeb se ef Oe KMUeCKNU SPL lef lh AeoMP1 eNUog INGe IN? 
14776 DIMENSION UC) 

1488. IMAG= ° 

L389. IF tORESCVED J eGts TOLFZ) GC TO 323 

1.39C e lf COACSCUCST LF «TOLLE LZ) GO Te 8 

£331 « ALPHA = -UCS)/U( 2) 

135? « RBETA> A1 PRA 

lore RETURN 

14946 8 IMAG= 1 

1495.6 RETURN 

14%Ge 1) RT= UCO)F€?> — 4¥L1) *U03) 

1437.6 Ie RE} 4es She Fo 

Lsu3.e 3 

LAdse Cc THE GUAURATIC OCES NCCT INTERSECT THIS CONST@PAINT © 
14206 € 

La01le PO IMAG= 1 

1402 6 RE TUR 

14036 ¢ 

14046 c THERE 16 A OOURLE FOCT. 

14056 Cc 

1400. 3Oo ALPHA = UCAIS(ULCLI*Ze) 

14)7. RETAS Al PEA 

1498. HE TUPN 

le 3%. c 

14106 « TwC REAL FONCTS EXEST. 

leil. € 

l4lee 49 NIS= CS2kTCPT) 

lLalse ALPHA= (€-U(2)-DNISI70C2e*UTLY?) 

1414.6 BETS= (-L62)401S 707-9001) ) 

l4el5e TF CALERA ol fe HETA) RETURN 

l+lBe. SAVE= ALPREA 

1417. ALPRA = RETA 

14ia6 HETA= SAVE 

141%. 2 TURN 

1470.64 [NO 

leele SUBROUTINE CONCEK(CGCMIER -KGMIN eMP2 ) 

1422. i 

1424.6 ¢ THIS SUeROQUTINE FVALULATES THE CONSTGAINT FUNCTIONS 
lee Ge € AND FINOS Tht SMALLEST VALUE IN GMINe ZFLAG IS 30TH AN 
1475. & TRPUT AND WUTKUT PAR ANATE Ie Tf ZFLAGet Gel INITTALL Ye THEN 
1426.6 € JINUY TRE NIKRLINZA® CONSTRAINT IS FVALUATED. PFLAG 1S SFT 
1427.6 Cc EQUAL TC 7 TF AKY CCNSTRAINT ITS NONPOSTITIVE e OTH RWISE IT 
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~~ 


LacedSe 
1L4°%e 
Laide 
lwitle 
letve 
lwitde 
1434.6 
14te 
14.3%. 
1437.6 
laive 
lative 
1447.6 
l44]le 
14%7e 
14446 
L444. 
l44%56 
144%.¢ 
1447e 
14446 
14a de 
14506 
1465le 
1452 e 
14536 
14546 
165% 
l4uhe 
1457. 
isSBe 
1453-6 
l46Oe 
la6le 
1457 6 
1463 

LenGe 
Lenoe 
l+utre 
14676 
146He. 
146%.8 
1470e 
1471e 
147e-e« 
La7 te 
La 746 
1475. 
LGaloe 
147¢e 
Lake 
lw /'., 
Las le 
1481 6 
14tl~o 
La Se 
LT4EHG, 
14735» 
14st 6 
14576 


Cc 
¢ 


RFEFUAING Zt tae 


~~ 
ro 


2 


) 


aw 
VO 


4S? 


IMPLICIT REAL *8 (A-HeO-Z) 
INTE GER Fie ILePE A SS oki epF7FLAGeP oldtin LD 
INTE GER? ISTYPFE et AelLE el Ae LE sPUNSLO (EZ) 


1.ul 

Y,IC CAD 
REAL ACGE 00} e61 290 3 CENTER (COND -E OMAX 2 SU INF 
INTFGE ee JHE 359) DI CHAI S52) KINEAS( 1397) 6 FCBAS(1 392) 
COMNMERSZNI aT HOLOSEL IS KCTS VSZO10O )FACCI 3210) 
CUMMON (SL OSTSGFIL (C10) eBFLSCIOPSEICIO eFACILIIM ee ICs C 
CEOMNOCNZEL CSTPZDOI C199 1CVsSNE CG 10) oF 105519) oF 2069019) 
COMMGONSLOE LZOLCESC2) eXxCL S02) 
CUAAONSENOCNS/GI (356010) 262040001 DIAC 359V6+RACaNO) 
COMNENZINRNDX27 JH eC I GMA SK INFAS Ss IDPAS 
CUMMINS SO FLS BT aN eo JI oMFLAG PDL ee POeMOI CINE OekKFUN SK UAC 
COMNONZUIMS THe& eM eOD eK MU ec KNUsRLL so Ml 2eMI LO NM INO, IND 


wMIN= Lie 


7#LAG= 9 
DO 7O° 1+ 14M 
k= yHCT) 
CG= Fe0 
IF (KAKUeCEG.C )} GC TE 15 
D9 19 J= feKRU 
vK= St KWU 
GC= GGt GICLeI*XC IK) 
CONTINUE 
3G GG+t PACI ) 
KX CLK Y= CG 
le €t® s6C. M} Ge te -2e 
fF (GGeGEeGMIN) GC TC 2¢ 
GMIN= GG 
KGMIN= I 
MP2s 1 


CONTINUE 
NO 49 I- te 
k= fICMACT) 
GG= CeO 
IF (KMUCEGLO) GO FE 35S 
of <¢ JS FsEKRV 


GE= GO C2CT os Jd*XC J) 
CONTINUE 
GG= CC# BCT) 
FLICIK)= GE 
IF (GGeGFeGMIN) GO TC 40 
GMIN= CG 
KGNMITR 1 
MP ¢ =] 
CONTIANG! 
IF (FDeFGeO) RETLRA 
OOT= Ce 
OOTI= Ce 
GG= Oe 


IF CINEQ - 2) SO 60s 
IF (KNUCE0e0D) GO 1G 55 
OG S$@ J= Lek NY 

JIK= Jt KMU 


NUT= OGI4 FICING eJ)¥XC IK) 

GG= FFLCING)= XCIKRG)*#(OOT+ ELCING) D 
TF (KMU.LEC*2O) GO TN FO 

WO &7 d= JekMU 

GG= (G4 CITING +3 d*xXO9) 
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UPS S. 
Jmeuy D 


Jy 


BGwewuwenwwhrWd 
~™’ 


BDNAUS 


et Oe me mee ee pe eet pt pew et 
vVSISTIS 
ee et eeeeeeee 


WwW 
< 
. 


1549. 
1541.6 
1542 « 
tHdS6 
1944.6 
15456 


ANAALRAM AVM 
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63 


me 


ae 


SUBRUUTINE IS 


SUBRCUT INE 


—E—E———— 


Gt) Te 6 
Kis2= KMU+ INP 
IF (KMUCEQe0) GO TO ES 


BG 3 J= 1eKMU 

ONTS DOT4 F2CINP oJ) *xK( J) 

NOTL= HOVi+ DECI AP es d**05) 

GG= BF OECINF)+ DOTI-“ *€KP2)*(00T+ FOCINP)) 
60 TO a 

[LF (GGesteGMIN) SE TURK 

OMIK= Cr 

MP2Q= ¢ 

ROTUIN 

END 


FIN FVALUATES THE FUNCTICN F(X)= 
DIACUXCVL)) FCF 1 eX (NU) ) > 
QTAC(CXCNG)) SCF O*X (MUDt FB) > 


C3CMUb+ CL*XEMU) ¢ 
<HORLI+ CF¥xX (ML) + 


SUBRCUTING FINE oY) 


THIS SOUTINE EVALUATES THE FUNCTICN WHICH THE FNOPOINT 


CURRENTLY TPYING TO FIND A ROOT OF. TF MFLAG 


EQUALS 1 THEN TRE CAST FURCTICNAL COMPRISING F 15 ONE DF THE 
INF QUALITY CUCNSTFAIATS=<—- Th! INEQUALITY DFTERWMINED PY THE 
PAVAMETERS MPC AND Jie 
IMPLICIT REAL #8 (A-HeC-Z) 
INTEGER PD ePOLeP CO eSS eFR oe ZFLAG OM ehOeCOlsSLLePLA 
INTEFCEE=2 ISTYPE et AeL Ee TAc TE ePUNSL C20) -1C( 8949) 
FL AL ACGCCO) -C OECD) o CMI RD Pe COND eFEVMAX « SUMINE 
NTECERSS JE 795))¢0F1GMACGES2) eKINEAS(2392)-.IOR AS(1397) 
CEANENSNEWTZ HOD C11) eo XO1 IV e201 I) oe ACC( 3619) 
CEMMCRSELCET/SEFIL (10) eFFECLOCLEICIO VS EACLOI*Ce I CelLC 
CUAMENZ PL CST2AZO1 C1001 Se OP1Go1O) of 14G019) 0209019) 
CLANMENAL NCENSZGI (759010) eG204000 15) e FAC S59)» + 720400) 
COAMON/INOATZ NUECT II eMURC 1S) »eNUC1 1) eMUC1 9) 
CO4SMONSSCALS 2Te Ne JS eoMELAG ND Le Pe ED oeMPID¢ INE Co KFUNS K JAC 
COMMPOINAZDIMZ LTHeN 0% eID eo KMU oe KNUC BL I 0 2L2eMPEi oe NM eo INQ 2 IND 


Nn 


INTEGER PLT e ALI se FPL 2 CLC 
DIMENSION F(DD9).YCNND) 
FCOC)= 4. 

KFUN> KEUNt DD 

KMthl= KMUF 1 

TeH= 49 

TF (KMULFC.I) GO 19 22 
OO Pf T= 1e¢KMU 

1= €§hEQSEGS2Y GE Fa § 
IF CELeEGeING) GO TO 345 
~R= j= theo 


DOT= Jed 
IF (KNU.EGQ.C) GO TOC IS 
DO 10 J= 16KKU 


JK= Jt KMU 


DOT= OET+# FICL esd? *=Y( JK) 

FCTRI= -Y¥CT) *(DOT+ ELCL))+ EFICT) 
TF CKMU et 920) GU 1 20 

O59 18 J= 1eKMU 


FCIR)= + CIRD* CCL es *YOI) 


1546. 
Ineo. 
Inds0« 
todsl« 
UnS27< 
$953 
Looe 
in55 
1556 
les? 
loss 
las 
1540 
Iso! 
15oe 
Lites 
15’.4 
Lae 
lor 
lnh? 
lbode 
1ASbGe 
15706 
1s7le 
{ Fee 
In73e 
174.6 
ln?Se 
1A7% 
tie 
In7Se 
i57S. 
I17306« 
inal « 
los32e 
Livi tbe 
15446 
Ilndtise 
19866 
1547. 
1598. 
158%. 
inGde 
18%. 
tr rae | 
WIS. 
iwms4e 
Inge 
Li 9%s 
BSP E-o 
1>F%e 
lo4-s 
LoOOD. 
LmHol. 
loJe6« 
1903.6 
Ltr be 
16056 
TUG e 
lH d7e 


e@ove@t@eoeeeee ee ove 


CUNT INY 

lf foLeeFCet) GO Tt 7C 
in] iINe®e® MY 

’ a* 1 “¥yt,oo 

twerr ft aM 

1f CINE .o#Getd GC TC 24 
ie €F e361) GM WW 35 

fas f= tea 

1 De 

1F(KMU. Ged? GO TO 39 

19) “s > 1 o*&™UL 

mT M614 FOUCIMB .sdevCs) 
(Pode =¥OLP SCOOT SO CINMP) D+ BF2C IME) 
ii AMUe! Co) GU TC @9 

a 4 LsaM’yu 
FC URd= FLIP + DP CIMB,JIFYOI) 
‘ Ti 

tuQOs 4 
CUMNTINY 

Tr (MEL AG et T2012 d KF TURN 

fe CvePUett+ oO) GE IC 129 
D0T=- De 

DOT1= be 
NOT?Z= OO. 

TF CINEOD = 2) B06 BO e¢ 
iF €KNUeF Qe GO {G&G 8€ 
WO 4 4u 1s NU 
JK= J+ <¥YL 
DAT= OL FLCIND oS ETT UK) 
F(OD)= OFLECINGI- YCINC) *(D0T+ ELCING)) 
TE (KMUeFCCI) GO TC 126 
p> #7 J= 9ekKMU 
cr 


F(OnI= FC 
GY TO 12) 
TE (KNULECeS)? GO TE 95 
DO 93 J= 16KMU 
T= CCT+ FICINP ede CI) 
ODOT O'Fl+ YN2PCTIANCeIP*YO) 
IKN= TN? KMU 
FC DOI= BPF2CINFI4 DCTIi- Y¥CIKAD*COCT+ ESC INS)) 
IF (MFLAG.LTe1) FETURRKR 
TF (MFLeEG SO) RE TURN 
= 2? 


Ti) + LI CINQe J )¥Y¥(5) 


1F (MFLAG Y £250 1¢€€ si Fd 
{F (MOD.ECe i} ¢c EE 24a¢ 

AF 4KNUCEECs9) GO TO 132 

NG 139 J= 1eKNU 

JK = Je FVYU 

FC OO)= F (OOF G1 CIs e SP*YOCUK) 
FCOC}= F CECT*+ BAC ISI 

RETURN 

1€ (KMU.LFC.O} GO TC 1680 

vO 1465 J= I9KMU 

FCOC= FCED)+ G2lCsseSI*VO I) 
F(DO)= F (CO)+ ERC) 

KE TURN + 
F(DDd= YOUN) 

RETURN 

NOT= Det 


Of 189 J= 140D 
NCT= CET# ACC(AssI*V( 5) 


1608.6 FCSCI= NGT- AY 


Lo dve RETUPN 
161¢. END 
lelle SUMROUTINE CLRIV 
loml?. C 
lolte C THIS SUP GCUTINE CALCULATES TRE EXACT JACOSIAN DF THE 
lela. c SILINFAK FURNCTICN Cee INEC FY BSLCONS. THERE MAY O# MAY 
loloe € NOT sf ANCTREP FUNC TICNAt APFENDED WHICH WE APE ATTEMPT— 
lole Cc ING TC MAKE SINDIANG wITR AREWTON®S METEDD. AL SU. A CDL= 
leal7. c UMN MAY oF CMITTFOC IF TRE FROPOINT ALECUORITHA FECUIRES 
lobe € THE CCRFE!SFECNOING VARTARLE TC RE FLXFC. 
LolGe C 
1629.4 I4PLICIT &EAL*A (A-F,C-Z2) 
lovle INTEGER FO PDL ePO2 eSS oP ep ZFLAG Pe Me COL Wile PL] 
loece INTFEGES ® 2 ISTYPF eh ALE so LAs TE PUN €02992.1CC899) 
lo2tse PEAL AC 4000) CC FO0) eCMIR ce CORD ERPMAX .SUMING 
loc4de INTEGE Re? JE C559 eo DTENA(GSE?) SK INEASOI 2673e1DORAS(1 397) 
lLo2Se COAMENSNC WIZ HOVOeT LI oe XC LO ZC 1LIIFACCH3e1 )) 
Loeie COMMCNZEL CSTSFELCLOD ePESCLOVSELCIOSESCIO) CC erelCoLC 
liv2? « COMMONZELCET2ZZDL CL Ie1O UI CS LC) oF LC Ge1I Dd oF 7209919) 
TO Se COAMONSENCONSZGI (459010) eG204O0C eo 1S) oe FAC 359) 6:1 0400) 
LisP%e CUAMUNZTIWT XITZ NMUFOLI) eMURCLO) »NUC1 OO) eMUCLD) 
lo ide COAMCNSSCALS BT edo SISoMELAG MDL ePe FO oO de INE CoKFUN eK JAC 
losle CUMMINJIIMZ THN eMe OD o KML eKNRU CIM Le PLA eMPl oe NY eo INGe IN? 
loJ?. INTECER COLJJ 
losse RJAC= KRJACH DE * HD 
1634. {fF (NFLAG AEG ed) KJAC= KYAC- ND 
le 356 NU 10 T= Leh 
lod6e QO 10 J= 1,NC 
lint7e 10 H(EseJI= Ce 
Ler d%e KAUL= * 4041 
lose KNUL= KNU¢ 1 
lLO40.6 901= OC-1 
loble COL JJ= J 
lo4aee IF (KMU-FC 63 GO FC 55 
lowe DO SJ J=1eKMU 
1044.6 TF (SeFO oJ eANDe MELAG CGTe1) GN TO 45 
16456 JM= JeCOLIJ 
lo4be TF tKMU-EG*O) GC TO 35 
1647. Ike= C 
Loar. Co © I= teKMU 
194%. TF COENEQCEGs23 GE YO iS 
WM5. Ir CPTehLQeINQ) GU TC 28 
165). 1S tee T= YRC 
lnoee FCTR «JM)= O1C1 eds) 
165 46 IF CPeNEso JF GC Fe 36 
LH 54. COT= O«9 
16556 1F (KNUSEGC6O) GE Th 25 
1OSte« ot "a & 1eKAU 
loS7e Khe= KMU4+K 
LH dhe 29 ODCT= NCT+# FL CL eK eXC KE) 
1059-6 25 H( ths IM) = HOCIReuM)- COT ELCT) 
16606 cai: 'Oo 39 
1nH5Le 28 Tec 1 
loole 39 CORT INGE 
lost. 35 If (FUAsEGsf) GCC TE $e 
1054.6 09 40 T= KMUT-eCh 
looble Ives T— KMU 
’ lores IF CINEOsEQet)? GO YC 37 


1607.6 IF CIVESFCeINGP) GCC TC 48 


1HG6Be , 

1HO6G. ? CECIVMR se JI-— FACIMP, JDHXOT) 

lm7U0e M 

Ilo7le 

le72e 

ln7te 

lm74. 

lo7%Se 59 CONTIAU: 

le7G. SS 1F (KNU.FCeO) GO TO 141 

1677. 90 79 J KNU1,.00 

lo7&. TF (Jet GeoJS eANDe MFLAGeCTe1) GC TO 65 

1e79 JM= J-COLJSI 

lore JMU= J-KMU 

lors) 1F (KMU.ECeO) GO TO Tz 

ins? 

loa 

lus4e + -Ge Ae St 

LONG e f 3 GC TC 58 

ltr4hre - . 

1lHA7-6 XCT)#F 101s JMU) 

LOM. C D 6¢ 

LoHde 

16706 

levle 

1o% 7 e 

lodw3e 

LersGe 

LOHYoe 

1oOVGbe *EQe 1) IRC= 1 

lersTe — (Gi zeECsO) GO TO 141 

lLH956 1) BO I= 1leKNU 

1629 L+K Mb 

173¢ TKM oe CF eJJS ec AND eMFLAG CCT el) COLIJIJ= 1 

1701 IKM. OC eJJS eAND. MFLAC CETL) GN TC 8) 

2 ING fTeF Cel) GC TQ 73 

I 
I 


eIKP} GO TO 7E 
1kO 
AMtil+ lk 
TKM - CCLJUJ 
(KMU.5Ced) GO TC 77 
no 75 * leKMU 
7S HOT K eKNC P= HOIK eo KNC)=— FSC T pK) eXCK) 
77 WEIR eKROV= HOTK oe KNC)=— FACT) 
60 To #®y 
78 IRC= 1 
ec CONTINUE 


IF MFl AG=1 THE INEQUALITY SPFCIFIEC BY MPD AND JJ 
OCTERPMINES TRE COTF ROW CF THE JACOBTAR MATRIX He 


141 ETF (MFLAGeL Tel) FETURA 
IF (MOPD.FC.0) GO TO 145 
GO TO 1°70 
IF CINFO-2) LEOs1ECLIEO 
bO 155 J= 1+20C 
tf CseGTekhU) GO TC £57 
HO IWsJ)d= CICINGs J) 
IF Cishis INO). GE TO 185 
DOT = 969 
TF (KNULeF Ce) GO TEC 154 


NOV ESWN KROL DNAS UVRK OC EWOLP eu 
See e_,@aeaqaeeece eee Beseeae11ee@est @dcernas @ 


SSN SNS N NN SS VSS SNS NSN WSS 
NEMS IN RODD DS me ee ee me re mem OO 


ee ee 


1 
1 
1 
i 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
L 
1 
1 
1 
1 
t 
L 
1 
1 
1 
1 
1 
1 
1 
1 
i 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
i 
i 


NS NN SSN NN NNN SNS SNS NYS NSN SS 


Sn 


NNN NNN SNS NNN 


OO 152 K= 1,.KNU 
KB= KMUF K 
DOT= DOTFt FICINO eK) *X( KE) 
H(DO»sJ)= IDCs) - OCT- ELL KMU) 
6& TE £55 
JK= J-KEU 
H(DCeJ)= -FICING -UK) *XCING) 
CONTINUE 
GG TO 1890 
KB2= FAP+ KMU 
OO 164 J= 1200 
fF €.506GT. KMU? GE TO _t6S 
HOOD wed = DECINEs JII-— F2CINP oJ) *K( KHZ) 
GO FG 754 
if {JeN* eKES) GO 1G LES 
HOT= Deu 
TF (KMULEQeI) GU TO ILEF 
DU 167 K= 1eKMU 
DOT = OUT+ F2CINE kK) *x(K) 
H(NCes)= -COT-F2CINP) 
CONTINU: 
Gd FO 1°79 
If €MEFULAG cE.T et) RETURK 
IF (MPC ecGeDd) RETURN 
IF (MFLAG = 2) 16002302249 
FF €MPEsEDe —1) €6 TGQ 210 
NQ 200 J= KMULEEC 
JM= J- KML 
HODGsJ)I= Cl (ISN) 
RETURN 
DO 220 J= 1eKMU 
H(DEeJ)= C2tJJeJ} 
VE TURN 
H(ODesJIJ)= 129 
PFETURK 
Ot) 259 F= 3s8D 
H(DD.eT = ACC (2.1) 
RE TUN 
END 
SUDPOL TING NORMOYseS16N) 
IMPLICTt?T RFAL*e( A-KeC-Z2) 
DIMENSION YC10) 


(ry*vC1) 


DSFNT 1S A SUBRLUTERKRE WHICH IMPLEMENTS DESCENT ON 
Nf iim CF JHE FUNCTICN WECSE RCCT I[S BEING OF TER- 
MINFD FtY NEWTON*S METEFOD. 


SUGPOLTINE DSENT CALERA .ERCRM)Y 

IMPLICIT WEAL#&& (A= CZ) 

INTEGES PD eH I1 eo PDA oe GSeRF pZELAG CF eiDeCOl eRL1, LF 
COAMONZ IIMS THet eM eDD e KML eo KNUCMH LePLeeMPleNMe I NQe INP 
CUMMERSNE WIZ ACTCoLLI oe XO LUIeZ0190).ACC(3.10) 
DEMENSTION YCLOVSF(1I0) 

DO 5S [=1.hC 


ee 


1799-6 
LAID. 
le dle 
lw d?e. 
LsO0Se 
LHed4oe 
ladse 
Lsdhe 
1éJI7« 
L508. 
1H)I9 
Lalo 
L411 
lale 
1814 
slGe 
LalSe 
L&T he 
Lol7e 
lolae 
1419-6 
LARD. 
lHiAl.s 
1427 « 
La2aise 
Lo2a. 
la2Se 
13?G6-6 
lda2Te 
LA]Q]K 
1aAeDe 
14.3. 
LA3le 
lasee 
laste 
1534.6 
laste 
lithe 
18376 
laste 
Lo3%e 
L44364 
1341. 
LAaG2e 
lUH436 
LH44e6 
14456 
1B466 
147 


5 ¥CRI= ¥*CT)=- ALPRACZC I) 
CALL FINE «¥) 
CALL NCEM(F 6 S20eDLE) 
IF (£2 eLTePNGRM) RETURN 
ALDHKA= ALFRAS.. 
TF CALPHA.CT. 2969071) GO TO 1 
RETURN 
END 


OL CCMF AND SOLVE ARE RCRROWED FROM FORSYTHE *S THF SOLUTION 
CF LINAEF SYSTEMS CF EQUATICNS. THEY IMPLEMENT A GAUSSTAN ELIM- 
INATICN SCREFMtE WITF PARTIAL SP IVUITING TC PRODUCE A LU DECGMPOSTION 
OF THE MATRIX Ae *SCLVE® FO RFORMS TRE RACK SURGTITIONS NECESSARY 
TO SOLVE UL * X = Be KONET FINDS SIGN( OFT A De 


SUBROUTINT DECOMF(NNSASLL) 

IMPLICIT REAL #8( 4-H.C-2) 

COAMEN/INTZ IPSC 20). KOE T eKOUNTeC I SING 
DIMENSION ACLI¢1T SO) eULCIO oe LOPSSCALETN(10) 


INITIALI#Z: IfSesULe AND SCALES. 
ISING= ) 
N35 S €= eK 


ULCTeJ))) 14262 


CONT LNUE 
IF €RCWNRM) 244,42 
SCALF S(T)= 1./7kCWRRY 
60 To $s 
CALL SINE(1) 
ISING= 1 
SCALESC(I)= de 
Ss CONTINUS 


GAUSSTAN FLIMINATION wITH PARTIAL PIVOTING 


DAES(ULCIPsK)) *SCALES( IP) 
{SEZE=876) 11.11.10 


lF CF 


nn 

= 
K pee ee Oe 
he ROH SO 


Trtain 


eonnOIGinn 
x 
Hamam TAZ NS 


ell fh ge OOO te ek 


eV 6 2 a 01 @ Oo Oo 4 BO eh Ol) Fe Se De 2 ES Pe. Os 6 AO Oe Oo Be CO 2 Be Se 8 8 Oe 8 2. eee 8 


ce 


-— 
No? 


ia 
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1a, 


ne 


KE>= 16 54K) 
PIVUT= UL(K2 6K) 
KP1l= Kl 


OC 1€ L= KPLeN 
Te= €PSC1) 
'M= -UL( IP eK) /FIVOT 
Ut (1P eK) = <-ENM 
OG t6 Js KPLsN 
ULCTE et) = ULI Pes) FMRUL (KA, J) 
CGNTING: 
CUNTIANUE 
KHR= IRPSCtH) 
TF (UL (KP ND) 19019619 
CALL SING(e) 
ISIANG= 1 
"eC TURN 
END 


SURRCUTIND SOLVE (RNSUL oft eX) 

IMPLICIT © AL €&( A-F.O-7) 

COMMCN/SINTZ LES 20) eKDET.KCUNTSISIANG 
OIMERSTUN YELC19219)2H010) 6x10) 


N= AN 
NOL= N+ 1 
EP=> FPS?) 
X(1)= ECITP) 
oO 2 I= 7.8 
IP= 19S(1) 
IMi= I-1 
SUM= 1249 
CO tL s= LelMtl 


SUs= SUM+# ULCIP es) *xCU) 
XCPI= EC TP )=- SuM 


IM= FRESEN) 
XON)= XONI/SUL CIP oN) 
IF CUL CIS 6NK)-GEeD) CD TC 1¢ 
KDET= -KDET 
, @ TRACK= 2,N 
I= NP1- It ACK 
I GCES (R-l)ecec sl 
TP= FPEtCE) 
[P1= I-41 
SUM= CeO 
CO 3 J= IPishK 
SLM= SUM4 ULOCTR es) *x0 J) 
CIV= ULCIP.slI ) 
{F (OIVeGE.9) GU TC 4 
KNET= -KDET 
ACT )= CXCL)-SUM)/JOIV 
KE TURN 
“ND 
SUMSOUTING FINGC IwkY) 
FO'’MATOCLX e* MATRIX wITH ZESC ROW IN OFCOMPNS=.") 
FLUMATCIXe "SINGULAR MATE IX IN DECOMPESEe 2°20 DIVIDE 
IF CIWhRY= 1) Lel o2 
MPFTTE (helt) 


IN SOLVE *) 
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CeWG Cw UDE) 
AL SH (A-+,0~-2) 


%> sPN1+P TP eSS» 
JHC 3590).N1GS 


SP ee ZFLAGeRSe Pe OMe dD! CHL I sAL? 
aqcs 

IS TYPE el ALE eIAe 
cu) 
Mie 


r 
ZV eKINFAS(1 302). IORAS(1 39?) 
TF ePUN LOO 290),I1C(83C9) 
C1STON {ac 
TSO eC BCID Del COND sERMAX » SUMINE 
CUMMCN SLW ,OPROD OY » DE SOP 080 3509 2XC 3593 .¥C3 
PAek ei MTS e COND oF RNA ee SUYVINE © LCNAM( 130202) NAME 
NTEMPE CP) eKINP STE TIM¢CJTIMSITINVeIJTINV oMSTAT OF Je T2044 IVINS TVOUTs 
LTOCNT se IMVEROSITSLIM, TIFFS Ze JCCLPeNRCW NCO a NF eNLELEMANLI TA, 
INGELEMe NINE eNUEL FM eNUETASNNFGNIeNLINE Se ISTY' 
SLACLIIO“’) eLEC2ZCOZ ).PUNT YE), 
HIPUNT eN CEC be SN QUAL NEP Its IF EAS*s LECRSH 
CC4MON LIVCHeLICHAsTFPIWT. IENFGeKOUTS 
CLUMMCN FAL4CC3)e TECBOCS) 
COMNMCN/LP1ZPT 0130: 
CCOAMONZ ALS (} 
COMMOERZ ILLS (1 
COMVYONZL NCONS “G1 €3! 
COMMOCN/INGKIZ NUE 
COMMEANSINEX?4 JHe? INBAS* IDR 
COMRCRZ> CALS SRP esd eM LAG eDO1Le FoF id eMPig TNE FUNeK JAC 
COAMUNSOIMZ M eID oe KML eKAUCRL Te FLAP eMO1L ENVY {Q. INP 
I= (wOD>-7) 
id CAEL CEFT LAC 
WOTTE( OE. 20) E 
22 FUGRMAT (1X68 é > LFFT £S NAW '.F9e5e* 
PE TLYEN 
3) WRITE (549) 
WAITE €2 041) 1 
WRITE (6.49) = 
fF (€KNU.EC.O) GO TF 5 
WRITE €e6e4?) OCG Tes) 
Sf fF €KPULEL6O) RE TURN 
WEITE C5047) (CG2CT ed) eo 1T= TeN)0d= 
RE TURN 
4) FORMAT 
42 FURMAT 
41 FORMAT 
49 FORMAT 
43 FORMAT 
6C TF (wb: 
Ga WRITE ¢ 
WRITE ¢ 
5st OOMAT 
= fe 


7c 


¢ 
e 
9 
i 
A 


~ Le 
Ie) 


Rg eace Rok a ae 


1968. 
1969. 
1970.6 
1971. 
1972. 
1I736 
1474. 
lv7Se 
L?Ge 
1977. 
1978. 
1Y79.6 
1980.6 
1S81l. 
19H 6 
l4Yedse 
19BS4. 
lw356 
198466 
1G87. 
1LU8B8.6 
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III. The HRA (Homotopy Retraction Algorithm) Code for Solving Equilibrium 


Problems 

III.1. Revisions to the Original Code 

The homotopy retraction algorithm is described in detail in 
Chapter IV of the report [2]. However, since the completion of that work 
the author has changed the HRA code to improve the domain of convergence 
of the code with respect to the user supplied initial utility levels. 
With these changes the code successfully finds equilibrium points for 
problems and starting points which resulted in failure for the previous 
version of the HRA. Also, solution times have been as fast, or faster, 
than before. This version of the code has not been tested on a sufficient 
number of test problems but the advantages over the old code seem to be 
so great that we should describe the code in its present form. 

By Lemma II.2.4, solving the equilibrium problem is equivalent 


to solving for 
if : a 24 
wo =: (Ors Ay 0 Ci Sim), P3 Ss ty 2 (Som), p) 
which satisfies 


f(w) = 0, weoD, (5) 
where 
n(w)we - A, (w) (t, (w) + v,) 


f£(w) = 


(ww - Cw) (t (w) + v,) 
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(1(w) refers to the first m components of w, etc., and 


D = {wly? tote, i€n 
pt zi +i pe =o = wi 

i=1 i=l 

mBr-d, y- o =0, i€n 
ft &--p =f 
Ts =0 

ce «0, i€nm 

w>0 } 


The original homotopy which is used to define a path to the 


solution of (1.2) is 
0 
F(w, 6) = 6f(w) - (1 - 8) (Aw) - A) 
where i” is determined by the solution to the auxilliary linear program 


(2). As part of the modification to HRA we changed the deformation above 


to 


Med 
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: 


Fiv, 0) = $86) - G- Ota) < 2°) 


Here t(w) refers to the slack variables of the first m rows, where 

m is the number of consumers in the economy. By the properties of the 
auxiliary linear program, it will always be true that en the value of 
t(w) at the optimum, is identically zero. This is because one can never 
achieve a higher level of exvorts p by allowing consumers to have a 


greater utility level than v Thus, we can simplify the definition of 


ic 
the homotopy to F:Dx[{0,1] > oe 


F(w, 8) = 6f(w) - (1 - 8)t(w) 


A convergence theorem using the path defined by (6) can be proved 
even though we will not do so here. 

Other amendments to the algorithm include a procedure for alter- 
ing the initial utility levels Vio i=i1, ..., m, if the processing of 
the auxiliary linear program indicates that the solution is unsatisfactory. 
There are three factors which can cause the program to change the initial 


utility levels: 


If the l.p. is infeasible, then v5 is too large. They will all 


be reduced. 


If some x” = 0 initially, it indicates that v, can be increased 


i 2 


without changing the optimal value of the objective function. For 


each i for which this is true we increase v, until Ay ¥ Os 


i 


3. If £,(w) <0 for some i, the initial boundary conditions are 
not satisfied. By reducing v5 for such indices i, we can 


increase f, until it is positive. 


i 


When we say to increase Vy or decrease Vv,» we mean that the 


following procedure is followed: 


a) After each solution of the auxiliary LP, a scalar TI 
is incremented by one starting with TI = 5. 
b) If any of the three conditions above occurs, to reduce 


Vv, we let 


sd v,*(1 =" RFEEY 


to increase ‘5 we let 


¥, =V¥04 1/TI) 


c) When TI reaches 20 and the solution to the LP is still 


unsatisfactory, the program terminates. 


III.2. Input Requirements 


The form of the input data for the HRA is exactly identical 


to the form for the BCA. See Section I1.1 for a description. 


III.3. Main Program 


The main program consists of 390 source statements which per- 
form many of the same functions as the main program for the BCA. The 
parameters are initialized and possibly altered by reading the name- 
list PARM1. The subroutines of LPM1 are called to input and solve 
the auxiliary linear program. Tests 1 and 2 of Section 1 are performed. 
The data is read so the C matrix can be constructed. Then the budget 
surpluses f. are calculated and test 3 is performed. The rest of the 
main program is essentially identical to that of the BCA except that 
different decisions are made after the ENDPNT subroutine has been 
called to determine which is the next cell (see Figure IV.3.1, Elken 


{[1977b]). 


Restrictions Relevant to the Use of HRA 


1. The number of consumers (IH) must be less than or equal to 10. 
2. The matrix D must not have more than 350 rows or 400 columns or 


more than 4000 non-zero elements. 


Of course, if one must solve a problem larger than these dimen- 
sions allow, more core must be allocated and the dimensions of the 
appropriate variables must be changed in every subroutine in which 
they appear. 


Before describing the other subroutines of HRA we will define 


the variables which differ from those in BCA as described in Section II.2. 


LMI. ELM a 


Variable Glossary 


In HRA we refer to one more variable from the blank common of 


LPMI. 


'N' if the auxiliary linear program 
is infeasible 
MSTAT = 
'O' aif the current solution is optimal 


‘'U' if the LP is unbounded. 


We use MSTAT to check if the current auxiliary linear program 
is infeasible. 

The common blocks Lyi, BLCST, BLCST2, LNCONS, INDX1, and INDX2 
are identical in both HRA and BCA. COMMON/INT/ contains the variable 
KEND, the number of times the ENDPNT subroutine is called. 

We include this variable to pass its value from the main pro- 
gram to ENDPNT so that we can test to see whether this is the first 
call to ENDPNT or not. In addition to those variables described in 


Section II.2, COMMON/SCAL/ contains the variable 


ERP = in 1 the number of consumers plus one. This is the dimen- 
sion of the consumer space plus one for the homotopy 


parameter. 


COMMON/DIM/ 
1, if the Newton tail routine is in effect, | 


ITAIL = 
0, otherwise 
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The description of the algorithm in Chapter IV of Elken [2], 
shows when ITAIL is set equal to l. 


The rest of the common blocks are identical to those in BCA. 


TIIl.4. Subroutines of HRA 
We will describe only the differences between the subroutines 


of HRA and those of BCA. 


1) Subroutine BLCONS: Calculates the coefficients of the bilinear 
functionals the budget surpluses so that f(w), (6), can be eval- 
uated in terms of the superbasic variables (see VI (4, 9) of 
Elken [2]. This subroutine is substantially the same as BLCONS 
in BCA except that in this case all the bilinear functionals are 
computed rather than the first DD, the number of bilinear con- 
straints which are currently satisfied. 

Subroutines called: UPAKC 

2) SUBROUTINE FINDP (PD1, IS, P): same as in BCA. See Section II.3. 

3) SUBROUTINE PIVOT (SS, RR), SUBROUTINE UPAKC (II), SUBROUTINE 
BSCNG (S, R), SUBROUTINE SUPERB (KEY, PK1, IS, (PD2, JS), and 
SUBROUTINE RECALC: same as in BCA, see Section I1.3. 


SUBROUTINE ENDPNT (JS, PD1, IS, NET): This routine implements 


the path-following algorithm described in Section IV.2 of [2]. The 


reader should consult that work to be able to understand this 
subroutine. Here we will describe some details which were not 


mentioned in Section II.3. 


The independent variables (superbasics) are placed in a vector 
X of dimension IH + 1 =IHP1l. The extra dimension is for the homo- 
topy parameter, 9. Theta always occupies the last position in X, 
X(IHP1). 

From computational experience, we know that the curve F+(0) 
is highly nonlinear in the last variable when the algorithm begins. 
That is, for the first tangential approximation in the first cell, 

8 increases very quickly, initially, but quite slowly thereafter. 
Remember, 6 increases from 0 to 1 as the algorithm progresses, although 
some decreases are possible. To remedy the poor guess at 6 which 
would result from following the initial tangent all the way to the 
opposite boundary of the cell, on the first iteration we move only 
one quarter of the distance to the opposite boundary. This rather 

ad hoc procedure seems to work quite well. After this initial step, 
the algorithm described in the work cited is followed precisely. The 
first tangent of the first cell is characterized by the fact that 
KEND = 1 (KEND is the number of cells traversed) and INFL = 0 (INFL 
is set equal to 1 after the first steplength is reduced). 

In this program, the integer variable MFLAG has a different 


meaning than in the BCA code 


0 when the tangent to the curve is being computed, 
when a "hyperplane subproblem" is being solved, 
2 when an endpoint in the opposite facet is being 
MFLAG = solved for, and 
3. when the equilibrium point appears to be in the 
current cell. Newton's method will be implemented 


as a tail routine on F(X) = 0. 
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SUBROUTINES CALLED: 


QUADS, DSENT, CONCHK, GTN, FTN, DERIVG, DERIV, NORM, DECOMP, SOLVE. 


Subroutines QUADS (U, IMAG, ALPHA, BETA) and DSENT (ALPHA, PNORM) 


are the same as in BCA. 


SUBROUTINE GIN (IFL, Q, U, F) evaluates the homotopy F(w, 6). 
When one writes this function in terms of the superbasic variables 


one has 


X(THP1) 


THETA’ f£ (X) 


KNU 
) G1(MUH(S), I)*X(KMU + I) + BA(MUH(J)), J = 1,...,KMU 
I=1 

-(1 - THETA)* 


The values of the variables of t indexed by MU (the first 
KMU in the vector above) are saved in the vector BLAM(10). These 
values are used in the subroutine DERIVG to be described later. 
The last component of Q contains the value of a functional 
determined by the value of IFL (MFLAG in subroutine ENDPNT). If 
IFL = 1, then 
IH 


Q(IHP1) = J} ACC(Z, I)*(X(I) - ACC(1, I)) . 
I=1 


If IFL = 2, then 


KNU 
Y Gl(JJ, I)*X(KMU +1) + BA(JJ) , if MPD= 1 


Q(IHP1) = 


} G2(JJ, 1)*X(I) + BB(JJ) , 
I=1 


The subroutine FIN is called to evaluate f(X). 


SUBROUTINE DERIVG (IFL, U, G, F): calculates the jacobian of the 


function evaluated by GIN. The jacobian is of the form 


0 Gl 


0 


DG(X) -|truzrsoco - (1. - mien ( 


SUBROUTINE CONCHK (GMIN, GKMIN, MP2): this subroutine evaluates 

the basic variables in both the primal and dual systems, saves the 
value of the minimum in GMIN and points to which variable it is with 
(KGMIN, MP2). From the theory of the homotopy retraction algorithm, 

it is known that the first IH primal and dual variables cannot become 
negative so we ignore these variables when searching for the minimum. 
Also, the primal and dual variables corresponding to the objective 


function are of no importance to the algorithm, so we do not compare 


these values when searching for the minimum. 


SUBROUTINE FIN(F, Y): evaluates the bilinear budget surpluses and 


stores them in F. This subroutine is identical to that in the BCA 
code except that no additional functionals are evaluated. That task 
is performed by GTN described above. 

SUBROUTINE DERIV: evaluates the jacobian of the function described 
in FIN. This is identical to the subroutine by the same name in the 
BCA code except, again, it is simpler because no rows need be appended 
because of an additional functional. 

SUBROUTINES NORM (Y, S1, NY, DECOMP (NN, A, UL), SOLVE (NN, UL, B, X), 
SING (IWHY), DEBUG (MODE): are identical to the subroutines with the 


same names in BCA. 


III.5. Sample Problems 

Below we give the output of the HRA program for the same two 
problems (MAS-Colell and Whisman) presented in Section II.4. The 
form of the input is identical for HRA and BCA so we will not repeat 
the input decks here. 

The output is very similar. The reader will note that the 
subproblem solved by ENDPNT is (1H + 1)-dimensional with HRA while 
the dimension grew from 1 to IH for the BCA. On these two problems, 
though, the HRA code required fewer calls to ENDPNT before equilibrium 


is reached. This behavior is typical with larger problems as well. 
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III.6. HRA Source Listing 


Below is a listing of HRA without the subroutines which com- 


prise the linear programming code LPM1. 


ETT NBN Ee I 


At 6g SENS AT 


erie 


oy 


OLAV AONE WIR ODEND DPW 
Qe 


ND ee et tte a ee ee 
eecret®t®eorettateereoetere 


Na 
ee 


ia) 


- 


AIVANIVAIVIANAAN 


REAL tA BP atl ns al 
INTEGER PDePO1LePD2eSeR 

INTE GSE*2 JHO35C) eC IGM 
INT® GES K2 TSTYPF sLA els 
DOUSLF PRETCISTON E( 890 
REAL ACA2OT) C(O SICD OS 


Ss 
* 
4 
’ 
Q 
M 
COMMON DSUMs PROD OY. DOE de YTEMP( 350), 
VA eS oe CME Ne SOND c0ERMA KS SUM 20). 
PNTEMP(2C eK INGO eI TIMeST I! JeTROWP eT VIN» I VOUT. 
ST TONTSTNVEE Ts ITRLEM SEE =M_ NE 
eT 350) 

3 

w 

w 


SNGULIMeNINFE sNUEL UM eo NU 


1 

e JCOLF eNRCW 9 NC: TAoNL ELEM eNLETAs 
“3 NNCGO J eNLINE So 
EL A(T 40? Veh ( 2002) oPUNG 
ELEUNC sMHEGT eNOUAL »NIPL F 
COMMON ITCHe TT CHA, IFPI: I 
SOMMOM TAC 4000) TE (999 


EAS» ITFCRSH 
FNEGsKOUT# 


DATA QN/IN Le 

COMMONZLP L/P C13 

COMMONS RL CST/ SFI 
1 
1 


: 
heal 
mn 


ef wOHe Tee OK 


COMMONZELCST?I 7D 
COMMONSUNCONE /G 
COMMOM/S INDX1LZ NU 
TOMMONSINDX2/7 JHe 
COMMOINSSCALZ STeNB 
COMMONSDIMS THeNoM 
COMMONSINTZ TPE( 30 
COMMIN/STILESZY TOLFZ 
SNIMENESTON V3(10) oF ( 
WAMEL TS T/PARMIZTOLFZ 6 VeTCNTRLs LECHOs 
TS TOMXsSTORD TCE CHCeolTHeleIP KOUTEsITLIME 


HRA 3 A CONE WHICH [MPL ZMENTS TFE FOMOTOPY RETRACTION 


we me 


) 
c 
1) 
Le 
1 
I 
. 
e 
) 


mee AOLOVLCe we 

OAKRZU ewe eo x 

2 4~ SUC: re a 

One KROPPFACNKO 

Ne HAN: GBR UAON 

oOnr-CCZe porervwr 
oH Ze OCe wm 
4aAnN56e Bruns ms 
Ime UPR OTO 
Teme NO~ww 


CLOF ZU“H4~Fe 
<e eo CH Z00K~ 


ALGIRI THM FOR SOLVING ECONCMIC EQUILI@RTIUM PROBLEMS « 


AUTHOR ¢ THOMAS Fe ELKEN 
SYSTEMS OOTIMIZATIUN LABCRATCRY 
PCOERATIONS RESEARCH DEPARTMENT 
STANFORD UNIVERSITY 


FOR DESTRIPTICN OF TRE ALGORTTHM AND DOCUMENTATION 
DF THIS PROGRAM SEE SOL TECHNICAL REPORT 77-AG. 


TOL ADs 1—C4 
TILCV= 10-C5 
TOL FE Zs 191% 
STPMXS 100 
STPLSHN= O05 
ITLIMF= 10 
TONTRL= 1 
TECHAs 0 
KOUTE= C¢ 


WITTE (6,PAIMI) 


97 


-* 
556 MUB= NFO 


Sle M= NECW 


: Sh. WRITS (6021) h 

57. 21 FORMATOS////, ¢ STRUCTURAL COLUMNS OF MATRIX. *) 

SA, 22 L3= MURF1 

5S. MUS= LA 14 

Ce K= 6 

le 99 25 J= LeeMUR y 

O26 IF (JeGTeNCOL) GO TO 24 

€3e K= Ke ! 

O46 CALL UNPACK(J) 

65.6 D0 24 f= 1yNROw 

Fe GLl¢leKdy= Y¥CT) 

A7 24 CONT INUS® 

ac, 23 CONTINUE 

Be 26 wRIT® (6027) LHeMUR 

7h 27 FIOMAT (47/77 4" TOLUMNS %o 140% THRUOGH *',14) 

Tle 90 25 LeleNSOw 

726 WITTE (692+) (GIL Tod) oJ=l eK?) 

726 2P FORMAT (1X elSFHe3 

The 23 CONTINUE 

7H, TF (MUR el Te NCUOL) GO TC 22 

re, 2° CALL NOPMAL 

7e CALL UNFAVL (0) 

Fe 1 60 71919 

Pea 39 sTOP 

thas 4&0 CALL INPUT (TACT) 

77.4 CALL OF BUG(1) 

7h 2 MFILLAG= C¢ 

7A, 3 1Toy= ¢ 

7B od T= Se 

72,5 343 IF CITFY ef Fe 10) GI TO 344 

7EAe WRITE (64 302) 

74.7 342 FOSMAT(* MUST TERMINATES TOU MANY TRIES AT INITIALIZATION 

79,9 1".7.* CF Vv WITHOUT SUCCESS.) 

FPS 344 IT?yY= JTRY + 1 

Pe ‘ CALL NORMAL 

7941 CALL UNPAVL (MFLAG) 

70,4? MFLAG= 7 

70,3 ¥ra FYUSe Le 

7%04 IF (“STAT .NFe ON) GO TZ 348 

7403 09 345 Ls telIH | 

7IeF 345 PCTI= (let -Le/TID * BCID 

7% 7 GO TO 243 

7ieP 143 90 3249 T= told | 

7969 IF (YCT) e6Te TAL3IO) GO TO 349 

Re (T= Cle® © LeATI)D * ECT) 

Roel WFLAG= 1 

H% —? 249 CONTINUE 
j 80.3 IF (MFLAG e3Te 0) GO TG 343 
; Ra, ¢ | 
; Boe CALL DE PUG(1) 
' B7. CALL SHIFTP (464) } | 
f Qn, KEND= ¢ 

2G, KENDSV= 9 

Ce IHt= T4e t | 


RAST Nas A Sa En aS LL RECN tt CE Bc nl neecemmanniesnets , aw a mates ernres eee 


926 N= NCO NROW 
936 NM= N@ M 4 
946 Lis Me t 
9%, MM1= Me 1 
Gh, LTSINVS 1 
a7. TE CLP 2000009) G9 TO 42 
~ 9%. READ CT el 7D) CSR CTV Tele Ih) 
\ : GG, 43 tc we? 
‘ ‘a 19 t3avec 
¥ 1O0le LCc(ld= 1 
; 102.6 70 AB K= lel 
— IC 3. KPYs <4 i 
| Ma, as RE ALC KINO S1OCLO) TSHae Lhe KK 
it 1SS 6 TF CISHReFQe1) GO TO €9 
10% {F (TS e8Ae-1) G) TC 42 
C7 IF (LL ei Qet) GO TI SE 
= 10% I™ (KK eNZ0O) GE TID 47 
10¢, [¢ Ques 1¢COu¢ 1 
; {tC. I1c¢Clcou)= Le 
™ Tile cq tcour= ACL) 
1126 ao TO 4€4§ 
YT1%e 4s? Ny 6a {s Lhe KK 
) tia. tcuus Cou 1 
i t1%6 tctrcsyu)s 2 : 
a Lite c(icours @¢L) ; 
117.6 ar CONTINUE 
r Tle. GD %:) aS 
| 119.6 a2 IF (KkeN5ed) GO TO 4a | 
12C 1cOu= tcnue 1 | 
121. IC(fCOUP= LL i 
1226 PEAY (551030) CCTCOU) } 
’ 123% GO TO 4S ' 
TPA. 44 TSAv= ICQue } j 
] 128 « 1 46 T= LieKK | 
oe LP. 1c cout 1 
L276 Ic(tco9use ft 
1 - 126, ae CONTINUE 
} 12°. GEAVCS ot 339) (CCI) es JEITSAVeOTCOU) 
/ 130.6 GO TO 4¢ 
= 1316 5° TE (KK ONE eG) GO TD 43 
13?. Icnus Icoue 1 
. 1336 TCC §COuds= LL 
. 134, C(TCOU) = SRP CK) *E(LLD 
j 13%-6 GO T) a¢ 
. 4 3b aa On SS Tt ite KK 
' 127. Tou [cue 1 
- 13k. COLCIU)= SRO (KD BACT) 
if 1 3%¢ Ic( 1tu)= 
14¢, 55 CONT INU 
Be 141. GO TQ 456 
14?. bial LOCKE Y= To NUE! 
143, 44% CONTINUE 
144, 1° (LCL CHIe7Q.4) GC TO 469 


—— 
- 
S 
> 
. 


Ea? ¢ MUI= 0 
N2TTE CH e41)d 
147. EL PiVeMweat(s7//77_¢ STRLC TURAL CCLUMNS OF MATRIXe *) 
14h. 4£2 = MWUG4l 
14°, “Us Let IH ~~ 
16C. «2 9 
1S1. Oy ©£S5 Js LBeMUs 


en 


99 


I 
i 
I 
I 
i 


CONTINUE 
CONTINUE 
WRITE (6457) LSeMUS 
FIQIMAT (4/7/77 6" COLUMNS %e 14s" THRUOGH 
1) 74 T= Lene Ow 
HRITE (6573) CGIC Tes) eo J=le TH) 
FORMAT (1X slOF190046) 
TZONTINUE 
N9 79 I= lem 
TR= JHCT) 
XXCTR Y= XCT) 

¥TEMo(T) 


©) GO TC 90 


APROV= OCEXYTEMPCIR) 
DSUM= DSUM + DPRID 
CONTINUE 


CONTINUFR 
r 


©) Go TC 105 


TAL. 
182. 
PSS 6 
1P4. 
185-6 
1FR. 
1A 7. 
188. 
1489, 
190.6 
191. 
139?.6 
19%. 
194, 
195 


Kaw SRA OMMUA 
xI2RnN SONS 


3 Hl BILINEAR PHASE OF THE ALGURITHMse VARIABLES XXC 1) eevee 
xxO'A) - = SLACKSe AND FI Cl) eoeoesFI(M) ARE THE USUAL PI"*Se 
xK(Me1) XOM4N) AFE TRE KYSe PIOCM41 Deco esPI (MEN) ARE THE 
OVAL SL MORE INITIALIZATIONSe 


wINANN 


er rr ra 


R 
| 
E 212. 120 CONTINUE 
% 2136 KNU= 0 
iy 2146 eave 0 aie 
99 4 = 
31a. Re eee (hk TMASETL0G%e0) GC FO 428 . 
SITs CALL SUPER S(1e lel 390) 
| | hie Gi TO 43¢ 
) | 219% 425 CALL SUPERAC Le-leleCeC) 
; 220 420 CONTINUE 
et) TF (KMUs°Q.9) GO TO 435 
: 2276 DO 433 I= 1+eKMU ; 
2236 332 VSCI)= PIC MUCT)) : 
224 435 IF (KNUe72e3) GO TO 440 
225.6 DO 438 f= 1eKNU 
asi tk= [# KMU ; 
227. 436 VOCTKR DS XXCNUCTD) 
226 CALL SLCING 
| 2206 440 CALL FING e VSD 
ve TO= 1 
ese WRITe lea CE CL) Tete th) 
e = ° i x 
sai Ke MOMMAD TINGE ENISIR, GOnGeT SURPLUSES 2 6¢719F13<6) 
re D5? 99 1230 15 LelH ae 
GAs If €°¢ 13 0G%e9) GO fo 139 
aE 257.6 ACL)= GOL) + BCLDZPLOL) 
« © 
‘ MFL AG= 1 
ai seas 140 CONTINUE 
361. IF (MFLAG «GT» 0) GN TO 2423 
| Saas ZFLAG= © 
= 6S 30 % aec 
DEK 4E 2 NTEMPQ1)= NT=MPC1)# UTI 
~ Dare “ FALL INVE@T 
{ 268 CALL PECALC 
v 2604 ITSTNVE 2 
7 a KEND= KEND+ 
ate COT TSInve [TSINV #1 ; 
~ Soe TE (KENG 2GT. ITLIMNS) STCP 
23 _ Lays ’ ry “t sis a eee 
| eae > FIAVAT / sia tlk Se, Ss a en na 6 bn OO ay 
is a7 5. bine oy: He Sofas? TH EMOPRT CALLE?) 
ane CALL FNDUNTOCSS PF OLelSeZFLAGeNET) 
. tee Ls ( (KENIDZ3)*3 eNEe KENC ) GO TO 468 
i SAC. CALL OF f8UG(4) 
a1 CALL DEAUG(T) 
“A ans. oP) SITWIT CiXet APTE® 'of3y6 NEWTON TTERATIONSS") 
As ¢ ee : M ~ 
224 > FLAGef£Qe2) GC TD 920 
i] daa: I (PD Ler ied cANDe f3eeGeNR}- 60. TC 920 
| a6. TF (O9Z) 30093100310 
Be Snes JOS RORMAT CIXC® THE "ofS. *Th DUAL VARIAGLE WENT TO 2880.°) 
be aa. 50 TO es 4) ts : 
2 = i . 
| Sate Sig FT JUMAT (IRL? THE Solte!TP CRIMAL VARTAGLE WENT TO ZERO. *) 
Bs 392% sBC IF (PD1e=Ge-1) GI 1) $20 
30a ¢ 2 cRO. THAT VARIABLE MUST GOTO 
| 7 2280 € NONGRET Cs AND TNE INCOMING VARIABLE 1S DETERMINED GY FINDING THE 
aa €  CARGEST ELEMENT EN Tee AOSROBATATE RCW OF Gis 
' 2966 € LARGEST ELEMEN : 
, al c 


298. PO= -1 


299% TRIWP= KINGASCTSD 
BAC e C4aLt. FINDP (Le TSeJCOCLP) 
BOL. TNNAMI= ICNAMCUCOLF et) 
Bierary INNAM2= ICNAM(JUCOLP. 2) 

q 3036 TONAM1= ICMAM(LISe1) 
304-6 TONAM2= LTONAMC( TS» 2) 
30Se WRITC(6 047°) INNAML » INN4SY2 6 LONAME » IONaM2 
BCA. SES TIOQMATC1Xe* VECIN:E %.2A4e! VECCUT: *.244) 
Rineary CALL PIVOT (UCILP.IS) 

: 3C8e TALL UNPAC K(JCOLP) 
3ODe CALL FTRANG(L) 
31Ce CALL WRETA 
Blle CALL 3BSCNG(JCOLP+IS) 
317. TALL SUPER T(Oe— Le JICOLP ele JCUP) 
3136 CALL PLCONS 
214. JS=2 £9 
Bld. IP CITSINV eGEe INVEEQ) GO TO 459 
31G6 sO TH. 462 
Jive Cc 
Bl ze Cc A OUAL VAXRTABLS wENT TC Z2ehGe TRAT VARIABLE MUST ENTER THE 
3196 Cc SASTS» ANS THE LFAVING VAFIAGLE WUST BE DETERMINED BY FINOING 
32Ce € THE LARGEST PIVOT FLEMENT IN THE IOBAS(IS)-ROW OF G2e 
321.6 Cc 
3A? G>2¢ Oe | 
3Pde JCOLPs 15 
32a, CALL FIND®2 Cele lSeP) 
225. TROWPS KINFAS(P) 
32. INNAM1= ICNAM( 1S 01) 
327.6 INNAM2= IONAM( 1S e¢2) 
BPR TLONAMIT= ICNAM(P,1) 
324. TONAM2= ICKAM(P,2) 
Che wWOTT® (64425) TNNAM1,INNAM2¢ITONAM1.IONAM2 
a21« CALL PI VITCJCOLFE .P) 
33BP oe ZALL UNP4&CK(JTOLP) 
33B3e6e CALL FTRAN(G1) 
334, TALL WRETA 
235e CALL ASCNG(JCOLP,P) 
3366 CALL SUPER ACO el «Porte ?) 
‘337 © ‘ £4. BLCONS 
BIR. JS= 16 
33%6 IF (ITSINVeSE eINVFFQ) GO TC 450 
34. GO TO 469 
341. S 
3426 ¢ ANOTHER 31LINEAR CONSTRAINT IS SATISFIEDe wE INCREASE IH 
3473.6 Lr AND WE HAVE 4A POINT IN wCIKDe 
34a. ¢ 
3456 220 WRITE (663930) KEND 
346.6 IAN CORMATC S/S," AFTEQ*414,* ENCPNT CALLS+e WE HAVE EQUILIBRIUMe®) 
347.6 WRITS (646220) KFUNsKJAC 
B4R. 329 FORMAT(Z." TOTALS?E SCALAR FUNCTION CALLS=*eIdeZe11Xs *JACOBIAN EVA 
349. ILVUATIONG = 61 407) 
35Ce CALL NE BUGO1) 
Ele 99 G35 Y= elIH 
3526 FBS VSCTI= MAC Te KOT) 
3536 WRITE (64940) (VSCT)eT=le ITH) 
354s WET TE (45e465) (PICT) e1Le to TH ) 
3556 POC FOSVAT (1X6* UTILITY LEVELS?2 "efF 1507) 
335A 745 © TVIMAT (1% e*OUAL MULTPLIERS: "eSF12e6) 
S57 6 7 CBS T= 1eN27NW 


102 


BSF Iv= JH(I) 
359. XCT)= XXC IV) 
3BoOCe ¥VOEIR= PICT) 
BAle 3695 CONTINUE 
3AP. CALL UNRAVL. (0) 
BE. TI. 
3A4as W219 FENGMAT C1E14a3 
JHSe LAE FORMAT CLIOFTS*4 } 
3666 ENO 
mS 3A7 6 SUPROUT INE A2t CONS 
3A. - 
JEGe 5 we ARS GUING TO CALCULATE THE JH ROWS OF COEFFICIENTS 
® 379. ¢€ FW THE B@ILINF AR FQUATIONS ANC TRE BGILINEAR INEQUALITY. 
3? ie Cc THE GASTC MATHEMATICAL STRUCTURE IS THE FOLLOWING: 
372? 6 e 
3736 c REL + OLSPTOMU)=— OLTASCEL (MU) )* Rice dit: + Fi) = 0 
3746 ¢ WP & D2*PITIMU)I=— DIAG X(NU)) CF 24%Pi(MU) + E2) = 0 
a 3756 Cc 
376.6 IMPLICIT FEAL*S (4-4,0=2) 
377. INTEGER PISPDLePO2sSeR pSSeRReZFLAGeRSeP 
37. INTE GERD JH0350) ePIGMA(GS2) sKINFAS(1392)- IDBAS(1302) 
379 © INTE GSE *D ISTYPE oLAck <0 lAy [2 ePUNSL Cl 20) ICCB 800) 
340. DQOUBL= PRECISION € (38009) 
3916 RFAL A(€40390).c( 80 C) eC MINeCCND sERNAX»SUM INF 
382. ¢ 
JAR, cOnNeTs CSUMs0PROIAD sOY 2 JE eC P08 (350) 0X (350 de ¥ (350 } e¥Y TEMP (350 Deo 
BAG es tA» Ee CVMIN, COND psERMAXs SUMINFE « ICNAM( 130202) sNAME(20) » 
3856 ENTEMP CAT) « KINPS I TIMeITIMe LTINVe ITINVeMSTATs LOB4e IROWPeIVINe IVOUTe 
3PF. ST TONT SZ INVES Os TTRLIMSIFFEZ « JCOLP Se NAC WeNCOL sNOELEMeNETAeNLECLEMsNLETAs 
BAT e SNGELEMs NINE sNUFLEM eNUET ASNNZGDJeNLINE Se LSTYPE( 350) » 
yaa, SLAC 1302) «L* (20025 ePFUN( 3) » 
BWAe ALPUNC Se eee ee IwWelFEAS,IFCRSH 
WWI. COMM ITCHe TTCHAsIFPLwls LENEGeKOUTS 
391 FaeMen TAC 40090) eI E( 8090) 
39276 CUVAMON/Z BLO ST/BFLC10) es SF2ACLIMZGE RC LOSE 2610) eC eo ICoilC 
3936 TOMMONZLP? ZET (1302) »XX( 1302) 
BGO, CAMMONZ AL OST2701610610)9202(9010) sF 149010) oF 219,10) 
WSe TOMMONZLNC ONS G1 0350 610) eG2(40041C 26 fA(350).RE( 400) 
JIA. OMMONSIND KLZ NUHC TO) »>MUH C10) eNUCIC) MUC1O) 
IGT COMPGNZINDXP7Z SHe OI GMAsKINEASs ICBAS 
BGRe TOMMONSSCALS BT e NB e JJ eo MFLAGeIHP1L Pa F Os MPD a KE UNS KJAC 
3Gve COMMONS DIM/S THeNeoMeKMU oe KNU o MP1 eNMe I TAIL 
400-6 Kiz= 9 
a2Cle K2= 6 
40°. TF (HM e& O29) GO TD 100 
a0 Se OG 8¢C *= 1 6«KMU 
404.6 K1l= MUCK) 
C5. fF (Kt1.GT.TH) 369 Ta 49 
aC6. IOD= KINGAS(KI 
: 407. TF (KNUsEQ69) GO TO 45 
40%. 99 49 T= LeKNU 
4CGe 49 Fi(Ke T= GICIDO1) 
q 41%. 6& F1(K)= BALIDOD) 
alle CALL UPAKI(«K1) 
4126 HO Ee = 19K™MU 
| 4136 SUM = °.0 
414.6 IF (M%UCT) el &e MD OSUM= YC(MUCT)) 
. AS 99 £9 Je 19M 
4166 IOJ= IOPAS(Y) 
A176 IF (ID) .FQ. 0) G TIO 50 


TES er 


4 


103 


A1% « 
8156 
427.6 
A?P1e 
42? 6 
427.6 
4246 
4?P5-6 
QP. 
APT. 
4226 
$296 
430. 
Sale 
A3?2Pe 
43 3 . 
4346 
S25. 
43. 
G27. 
43h e 
4BBGe 
44.4 
441.6 
442. 
A346 
444. 
B45. 
445.6 
447. 
ar%e 
44S, 
ASC 6 
S31. 
AS2e 
4536 
454. 
a5Se 
a5 6 
457 
452, 
45S. 
GEC 6 
&€ 1 e 
4E2e 
4636 
464. 
4556 
Abt. 
GET. 
4AM» 
445. 
&7C e 
471.6 
4726 
&T3 
47% 
4756 
a7E.s 
477. 


aAINnA 


AAANY 


wn 
5 


79 
at 


PICK) IS BAdICe 


10¢ 


a 


NN 
wQ 


Naw 


140 


145 


OSUM= DSUM+ YC J)*G2(IDJ ol) 
CONTINUE 
31(K er) = OSUM 
CONT INJVE 
OFUCKI=S Ce 
oO FO Je 1.™ 
IF (TOPAS(J)2F O09) GO TC 70 
AF LOK) = BF LIK) + YC J) *BSB(TDBAS(J)) 
CCNT INUF 
TONTINUE 


NOw THE S=SCOND TYPE OF EQUATION IS CALCULATED 3ECAUSE 


TF (KMNUeTQ29) RETUFN 
DN 146 K= 1e¢KNU 

K2= NUCK) 

TF(K2eG7T.1H) GO TO 145 
TOD= INEAS (K?2?) 


TF €KMUC 229) GO TE 123 

Dd. 120 T= lekKMU | 
F2(KeIv= G2CTDDe!) 

E2(K)= 8S( TDD) 
CALL UPAKC(K?) 
IF (KMU.FQ. C) GO TO 132 1 
AO 139 f= 1+KMU 
D2CKsT) = 029 . 
TE (MUCT) el Se M) D2(Kel)= YOMUCT)) | 
SF 126 J= tee 
IF( ISBASCJI)62Q00) GO TC 125 
N2(KeT)V= DA(KeTI+ YOIV*G2C TOBASO I) oT) 

CONTINUE 

CONTINUE 

SF2(KI= De 

797 145 J= Lem ! 
t= (IOPAS(I)e£Q20) GO TO 149 

AFQACK)= BF PCK)I+ YC J)*BECICEAS(J)) 

CONT TNUE 

IF (K2eNF.LTH) GO TC 145 

CONT INUF 

QETUEN 

=ND 


SURFOUTING FINDPC POLI SP) 


THIS SUBFEDUTINE CHCOSES Tre VARIABLE TO S8ECOME IMPLICITLY 


BASIC AS THE ONE WITH THE PIVOT ELEMENT LARGEST IN 
AGSULUTE VALUE 


TMOPLICTT KEAL*S9 (AH, 0-2) 

TNFEGER PD ePD1e PIDs S4ReSSe Re ZFLAG oR SP 

INTEGER *2 J4O 350) ,01TGMA(SS2) eK INBAS( 1302) eIDBAS( 1302) 
SIMMONZENTUONS/G1I 0350 619) 6G2(80001C de BA( 350 )e BBL 400 ) 

COMMONZINIXIZ NUH(C10) ¢MUH( 10) eNUC IC) »MUC1O) 

COMMONS INDX24 JHe DIGMA,KINEAS, IDBAS 

TOMMON/SDIMZS THeNoeMeKMUe KAU o MPL oNMolI TAI 


Be ff 

VE CHOP sttier th} Go Tt 39 
TOD= KINBAS(IS) 

99 20 I= 1 eKNU 


104 


- 


—_  —_—__— ———— iro 


a7Re IF (PeNE&E-9) GO TI 10 
479-6 AIG= DARS(GICTOIeI)) 
43C.e P= NUCT) 
aBle Ge F2) 2e 
ABD 1C COMP= PANS(GLCIDI,1I)) 
4383-6 IF (COMP.LLE.3IG) GO TO 20 
484, fStTG= Comp 
ARS. D= NUCIT) 
ABR. 2° CONTINUE 
GET YF (BIGeLTeTALFZ) P= 9 
4BP. PETURN 
A396 30 1NO= IDBAS( TS) 
Aga. DN 49 J= 3 .KMU 
AQ1e Fe Ce&ghe 6) 32 F238 
497% e PITS= M4lS5(G7EIOO.J)) | 
AG3e P= MUC 5) i 
4946 35 TOMP= DAdsS(G2C1IbDeJ)) : 
4956 TF (COMP.LEe8IG) GC TO 40 ; 
ADFe RIG= CiyMP 
457 6 P= Mud) 
4&S3e 40 CONTINUE i 
4996 IF (2IGeL TeTILFZ) P= D0 : 
509-6 RETURN | 
5Ole ND 
5027 « SUAROUTING PIVOT (S3exXR) 
5Cze € 
S04. Cc THIS YCUTINE PERFORMS A PIVOT CN THE PRIMAL SUPER= 
SSS6 S RAST COC COLUMNS IF POeEQele OUAL IF PUD ecEQe—1e THE PIVOT 
SCGe ¢€ BEINGS COLUMN SS INTO THE EASTIS ANC CCLUMN RK OUT OF THE 
S97 e e AASISe 
SORe Cc 
5C 3. IMPLICIT 8c alL*8 (4—HeI~Z) 
51C.e INTEGER PO ePDLeP IA ¢SeX esSeRreZFLAGeRSeP 
Stile INTEGFEX2 SHC 3BSO) -DIGMA(SS2) eK INPAS (1302 )¢ IDBAS(1302) 
Sie INTEGES *2Z TSTYPE SLA cLE sp TA, IE ePUNSLC(20) »fC(800) 
Sie IOUBLE PRECISION €(89C0) 
514.6 REAL A(4080) «CO 300) se CMIN»s COND cERMAX eSUMINE 
SiS. c 
SrA. TCOMMIIN DSUMs IPRODO sDVeNDE oOP 6:3635C) »X(35))+VY(350)}+YTEMP( 350). 
Sive LAs Fe OMI Ne COND «= RMA Xe SUMINFE oI CNAM( 120202) eNAME(20) © 
Slee OPNTEMP CAC) se KINPe IL TIMeJITIMe LTINVeITINVeMSTAT oe LOB se IROWP eI VINe IL VOUTe 
5ide BI TONTeINVERGs I TRLIMeIFFEZ se JCOLPsNRCWeNC IL eNELEMep NETAsNLELEMsNLETAs 
520. ANGEL EMs NINFeNUEL EMsNUETAsNNEGDJeNLINESs ISTYPE( 350) 6 
SPt 6 5LA(1 302) «LEC? 902) ,FUN(3) 6 
See FIPUNC eNDAGITeNOUAL eNIPILWeIFBASSIFCRSH 
SPS COMMON ITTCHe TT CHASIFPIWTs IFNEGeKOUTS 
Soh. COMMUN J4(6009)e15( 9900) 
7 5256 COMMONZL PLZPT (C1 3902)eXK013202) 
S266 COMMONS LNC ONS/G1 6350010) 2620400010) +3A(359 d+ BBL 400 ) 
S286 COMMONS INO XYZ NUHCLO) +> MUH(C 10) ¢eNUC1C) eMUC1O) 
52e. COMMONS TNDX24 JHeOLGM4sKINEASs ICBAS 
S306 COMMON/DIMZ LHeNeMeKMUsKNU o MPI oNMOeTTATIL 
7 S31. FPS= t€ .**(-1 3) 
S326 Pes KINGAS(RR) 
> Sea's IF (355 eGTe IH) GO TO § 
524. fF (NUM(SS)eNE el) GO TC EC 
Save S§ CALL UNPACK(SS) 
° 5366 CALL FTRAN(1) 
S37 YTEMPCRE)= -1e/YC(RS) 
be 538. 00 160 f= tom 


105 


rust sno et NS NE tS = ASRS SA NS RR 


DINAN 


- 


ela 
ee 2 
LO Oo 


129 


TF (LeNFekR) YTEMP(T)= VYOLI*SYTEMP CRA) 
CONT INUE 
60 TO 3°f 
JS= NUW(SS) 
YTEMP( RE D= 1Le/GI(REJS) 
90 ?S l= Lem 
IF (CTeNf eR3) YTEMP(ID= -GICTeJS)*VTEMPCRB) 
CONTINGU® 
TF (KNJ 05009) GO TO $5 
90 59 J= 1eKNU 
TE (NUC J) 66 O6SS) GC TI 56 
a 


de 
36) 
VeYTEMP( I) 
) 


V= ACRF) 

SA(PS3)= Co 

9 60 T= tem 

BACTI= GACI)+ v¥YTEMPC(TI) 
BAC FR)= -BACRB) 


we ARE NOW PIVOTING CN THE DUAL SYSTEMe IF THIS 
NOT AN IMPLICIT 3ASIC=TYPE PIVOT. wE CALL SUPERB 
CALCULA T" THE PIVCT COLLMNe 


IFL= 9 

RA= IDORASCSS) 

IF (RR «GTe IH) GO TO 7C 
IF (MUH(RR)eNE*O) GO TO BO 


YTEMP(RB)= 1e/G2( RBIS) 

90 85 T= 1eN 

IF (TeNEC?B) YTEMPC I= -G2¢f eJS)*YTEMP(RI) 
CONT INUF 

TF (KMUe7O09) GO TE 115 

99 1190 J= 1LeKMU 

TF (MUCJ)eFQLRR) GO TD 110 

= G2¢F 


G2(PdedJ 


ce 


wue 


.# V*VYTEMPCT) 
Bed) 


" 
f 
n~e 


w7ncCUC~ 
DIS tew DT 
a 


© 

leN 
43CT)= PACL)+ V*YTEMPCT) 
SACRPA)= -SACRB) 

20) RETURN 

FB 020060 e-d oRR) 
RETURN 
END 


SS _—_— o_o 


te bas ateabWAlaceh md ct se piers ae SA Se 


& ’ 


Blle 


ia) 


MAAN 


AD NDNAANNADAY 


PMPLTCET REAL*R (A-H.9-2) 

INTEGER*®2 JH( 350) sCIGMA(952) eK INBAS(1302)+1O0GAS (1302) 
INTEGER? I[STYPE,LAsL He LAs LE sPUNSLC( 20) 5101800) 

ONURAL® PreK CISTON E¢s9CQ0) 

QE AL 4(4700) eC 40C) eT MINs CUND sERMAXs SUMINE 

COMMON DSUMeDPROD DY ANE sOF 620350) oX(353 90 ¥ (350 do YTEMP( 350 De 
LA e& eCMIN + COND eERMAXs SUM INF e LCNAM( 1202 02) eNAME(20) » 
PNTEMPCZO0) oe KINPe LTIMs IT IMe IT INVe IT INVeMSTATs [OG4e IROWPe ITVINe LVOUTs 
SI TONTeINVERQs TTRLIMSIFFEZ ep JCOLPeNRCWeNCOL eo NELEMs 
&NGFLEMe NINFeNUEL EMe NUETAsNNEGDJSNLINESS IST YPC E( 350) » 
SLAC 13292) L&E (2092) -PUN( 3) » 

ETPUNCeNDEGT NOUAL eNIPI We lF 3A4Ss LECRSH 

COMMON ITCHe LT CHAsTFOLWTs LFNEGsKOUTE 


COMMON TAC SDIODV ee LF CL 3000? 
COMMUN BLCSTS BEL (10 e BFZCISISEICLOD E201 O) oCeICeLC 
CSCOMMUNJOIMZ THeNeMeKMU a KNU o MPL RAST TAIL 


vers Oe 


° 
o 

ra) 
o 
2 
4 
~ 
Zz 
c¢ 
_ 


o 
2s 
2 
Zz 

=~ 


WE TURN 
EN. 


THE FOLLOWING ROUTINE MAKES THE CHANGES IN THE INOEX 
SETS NECESS4=Y EVERY TIME A BEASTS CHANGE IS MADEe 


SURROUTINEGE ASCNG(Se2) 

INTFGFE* 2 JA(3250) eDIGWA(952) eKINGAS(1302)0e IOBASC13C2) 

COMMONS INIK2AZ JHe OI GMA cKINEASSICEAS 

INTE GEE SeSSsAVeReF SAV 

RSAV= KINGAAS(R) 

JHERSAV)= S 

KINSAS(CF)= 0 

KINBAS(OS)= 

SSAV= ITORAS 

DIGMACSSAV) 

INVP45(5)= 4% 

TOSAS(R = SSAV 

RETURN 

FND 

THIS SURROUTINE CAN OC THREE THINGS DETERMINED BY THE 

PARAMETER KEV, IT CAN ACO A COLUMNe REMOVE A CULUMNe OR 
BOTH ADD AND SEMOVE CCLUMAS FSOM TRE FRIM&aL CR DUAL SUPER= 
QASTC COLUMNS « DEPENDING ON WHETHER KEY IS Le 2e OR Oe RESP 
ECTIV =LYe POL TS 1 OR -=bL CEPENCING ON WHETHER A PRIMAL OK 
DUAL COLUMN I RE ING ALDEDe Fo2 IS A SIMILAR FLAG FOR THE 
C™LIMN BEING wnEMOVEDe IS AND JS INDICATE THE PARTICULAR 
COLUMN TO FE ADDED GR REMOCVED FROM THE GROUP OF COLUMNS 
SPECTEIEN BY PDL AND PD2e 


107 


NETA»NLELEMeNLETAs 


| 
/ 


a 
7 


459° 


64le 
66276 
6636 


644 
647. 
BFR 
FAD. 
ATC. 
S71. 
6726 
ATR. 
674.6 
BTS 6 
57He 
ATRe 
S736 
630-6 
6Ble 
FA2 6 
6B 3 
AAR. 
6OASe 
HBF e 
SBT 6 
6386 
689 
69C. 
691-6 
6926 
597s 
AI4. 
6956 
EGAe 
697 « 
BIR. 
6996 
7C06 
791. 
POP © 
7O36 
7046 
7056 
706-6 
797 6 
TOR. 
7C9e 
71C.e 
Tlie 
7126 
VFlae 
7i%.s 
T1556 
TIE. 
T17 


“ 


20 


25 
2¢€ 


SUBROUTINE SUPERA(KEYePOlseISePO2eJ2S) 
IMPLICIT REAL*S3 (AmHeO-2) ) 


INTEGER 2M ePIDLePD2o Sor eSS er eZFLAG eRSeP 

INTEGER*2 JH 350) eDIGWA(GS2Z) eK INHAS(1392)+ IDBAS(1302) . iJ 
INTEGER*® 2 LSTYPE eLAcLE olAs TE ePUNeLC (20) e1C(800) 

QOUSLE PRECISION E(8000) 

PEAL A(42000) 6€(300)eCMINe COND cERMAX oe SUMINE 

COMMON DSUMsDPRON DY e DE.CPeB(35C de X( 350 )2¥(350) »sYTEMP( 350), § 
Ase eCMT Ne COND 0 2PMAXe SUMINE oe ICNAM( 1302 02) eNAME(20)> 
ent’ EMP(2C)sKINPS ITIMsITIMeITINV pJTINVeMSTAT el OBJ ce TROWP eI VINe IVOUTe 
BYITONTeINVEF Os ITRLI Mel FF EZs JCOLPeNRCW eNCOL eNELEMe NE TAsNLELEMsNLETAs 
OaNGELEMe NINE eNUELEM sp NUS TAe NREGDJ eNLINE Sol STYPE( 350) « 
£LAC1302 deLE (2002) ePUNCS Do 
SLPUNC eNDF GI sNOUAL eNTIPI Welf EASs LFCKSH 

COMMON ITCHe ITCHAs TEP LaTe IFNEG eKOUTB 

TIMMON TA( 4000) eTE (32909) 

COMMON/LOLTZOTO L302) e XK( 1202) 

COMMON/LNC ONS /G1 0350010 )262(400010 )e BAC 350 ) + 68( 400) 

COMMONZINOXUTZ NUHO1O) eMLH(1LO) »NUC10) oMUC10) 

COMMUN INDX2/7 JHe C1 GMA eK INBAS»s LDEAS 

COMMON/SE CALS ABT e NBe tse MFLAGe IHF] oP oe POe MPCe KFUNS KJAC 

COMMON/DIMZ THeNoeMeKMU oKNU eo MPL op NMel TAIL 


EOS= 16 e** (15) 

{F CKEV eK et) GA TO 59d 
IF (P92 2.5T 00) GO TC 2) 
KM= KMU 

KMU= KMU~ 1 

TF (49S eGTe LH) GO TO 4 
MH= MUH(J5S) 

MUH(JS)= 0 

IF (MHeEQe KM) GO TC 6 
99 5 Y= MHeKMU 

IM= MUCT #1) 

MUCT)= I” 

MUHCIM )=  iphbshicaiie 1 

29 §& Sz le 

G2UJel)= é2tu. 182) 
CONTI NUE 

MUCKM)= 0 

99 7 J= 16N 

G2 SeKMI= O29 

IF (KEY e£Q0e2) RETURN 
69 TN S¢ 
KN= KMU 
KNU= KNU- 1 

IF (JS eGTe IH) GC TO 28 - 
NH= NUH(JS) 

NUH( JS) = 9 

TF (NHeEQ*KN) GO TO 25 
OO P® T= NHeKNU 

JN= NUCT 41) an 
NUCT)= JN 
NUHCJN) = NUHCUN) = 1 

90 25 J= LemM 

GICIoT)= Gl IeI41) 
TOMTT NUE a 
NUCKND= 9 


OEE RE 


NT eee a are 


7206 


7227.6 

ee 7236 
72a, 

| PPS 
} THe 
727. 

j T2B 6 

729 
7306 
P31. 
PFI3IZQ6 
TIAde 
Pike 
7380 
Tafbe 
7T3A7e 
7336 
T3BSe 
FaCe 
Tale 
742 6 
7436 
T4%e 
7AaSe 
7466 
Ta7. 
T4596 
7T4a%e 
75060 
75le 
TH? © 
7TS36 
7546 
7TSS86 
7TS66 
Ts? oe 
THR 6 
FEO, 
TOC. 
Tale 
TH26 
7636 
TEA» 
TESe 
The 
TH? 
THB e 
TE%e 
T7%C« 
771.6 
7T2e 
W736 
774.6 
VIE 6 
77h. 
TTT ¢ 
778. 
T7796 


nan 


GAN 


andan 


woMWWAN 


50 


SS 


490 


70 


73 


Is 


74 


9° 


10% 


TS PARTLY CL*EITNVIBD 
~CL¥®INVIA)D 


THE 


110 


IF (KEY e602) RETURN 
NOW Wo ADD THE SUPERBASIC VARIARLE SPECIFIED BY PDI sISe 


TF (POL-LT.0) GO TC 70 
KN= KNJ 
KNUS KNUEe 


IF ('SeG6Te TH) GO TO 5S 
NUCKNUD= LS 

NUHODS) = KNU 

CALL UNPACK(TS) 

CALL FTSAN(1) 

DQ &€f Ts tM 
GICIe<NU)= -YCD) 

RE TURN 


NOw *®£& ARE GRINGING 4 CUAL SUPERBASIC COLUMN INe 


KMU> KM¥UF I 

TF (€135¢e¢G6TetIh) GU TS 73 
MUCKMU)= TS 
MUHMOIG)= KMY 
IF (YSeLEeM) GO TO 110 

THF ENTERING VARIASLE 1S A ZETA. THe cNTER ING CULUMN 
& ROw OF INVCR) AND A BCWw CF INVOR)*De 


1s” 


c 

2¢ 

ONT INUF 
9100 v= 4M 
F (LO8ASOCYU 
ALL UNPACK 
IT= 9eC 

9 9C I= 1M 

IF (IDPASC12eEQ09) GO TA 99 
NNT= NOT# YOL)*xXOD) 

CONT INUE 

G2CIP 8 ASTI) +KMU)= DCT 
CONT TINUE 

SE TUPN 


a~t 


THE FNTERING VARIABLE IS A Pie THE ENTERING COLUMN 
AND PARTLY C24 Ci *INV(B)¥*D. 

TS 2ART OF THE BASES INVARSE CURRESPONDING TY 
X=-8AS1C COLUMNS AND TEE TeBASIC FOWSs 


10%= KINRASCIS)D 
30 170 T= Lem 
YC Co 

YCITIOM)= Ie 

Catt "TRAN 

CALL SHIFTR(I3-2) 


109 


PAC. 29 140 T= 1M 

TAL. TF (IOP ASCID«EQeO) GO TI 1490 ad 
7TA2 6 S2(IOBAS( I) eKMU)= X(T) 
7H36 1¢C CONTINUE a} 
784. 29 160 J= MPL eNM 
TAS. IF (IOBAS(J)2eEQe9) GOA TO 160 
TAG. CALL UNPACK(J) '? 
7e7. 297% 040 
TAR, AQ 150 T= tom a 
THQ. IF CIDRASCT)2£&Qe9) GO TC 150 
IGM, OOT=S ON7+ VOLIeXCTD) 
FOL. 156¢ CONTINUE - 
7926 S2CINDBAS( SD eKMUD= Y(IS)I* COT 
79360 160 CONT INUE éh 
794, KE TUEN 
71956 FND 
ao Lae SUSFOUTING RECALC 7 
POV € j 
74H © c THIS SUBROUTINE RECALCULATES THE SUPERSASIC COLUMNS a 
70S 6 C USING THE CURRENT BASIS IN GR FORMe It CAN ALSO ADO 
oats 3 A VARTARLE AND COLUMN TO THE SUPERBASICSe 
Ole = 
BO26. IMPLICIT REALBR (A-4-.0-2) 
AOBe INTEGER PD +PD1ePDO2eSeR »SSeRRoZFLAG eRS oP #) 
Be INTE GER*2 JH( 3530) eCIGMA(992) eK INEAS( 1302). IDBAS( 1302) 
IC 5e INTEGER 2 TSTYPF eLA ele eLAe IE ePUNeLC(29) ef C (A000) 
Ake OOURLE PeFECISTION £€ (8900) - 
ROT. MEAL A(4200) eC 0400) eC ¥INe COND pERMAX pe SUMINE 
BOk e * } 
PC. COMMON DSUMeDPROD AMY s DE se CFs B(350 ) eX (350) e Y( 350). VTEMP( 350)6 4 
B10. LAs = oe CMINe COND cE R MAX e SUMINE oe LCNAM( 1302 02) eo NAME(20) 0 
Alle PNTEMP(20 De KING ST TIMeITIMe ITINVedSTINVeMSTATe 1OBIe [ROMP eI VINe IL VOUT. . 
S126 AT TONTeINVERGs I TRLIMSI FF EZ oe JCOLPeNRC wWeNCUL eNELEMe NETAoNLELEMeNLETAS 
ALMe ONGEL EMA NINE eNUEL CM eNUE TACNNEGUJ eNLINE Se I STYPE (350) 6 
314.6 FLATL IO?) eL (20902) ePUN(B8 De 
A1lS. FITPUNCeNDEGIeN UAL »NIPI weIF BASeIFCRSH 
Ble COMMON ITCHe ITCHA, IFPIWwle IFNEGeKOUTH 7 
ALT. COMMON I[A( 4000) -1E( 3000) 
BIR. COMMON/L PL/SP 101302) 0eXKKC 1202) 
BP0e COMMONS LNCONS £61 450019 )0G2(400 01 C) - BA 350), BB 400) 4 
APL. COMMONZINOX1Z NUHO 10) eMURC IC) oeNLO IC) eMUCLO) 
A226 COMMONS IND X27 JHe CIGMAsKINEASe IOBAS - 
3236 TOMMONSSCALZ AT eNYe Ite MFLAG eIHPIL eP ePD eo MPO eo KFUNS KJAC 
S246 COMMONS OIMZ LHoHeNoeMeKMU eo KNU e MP1 eNMe ITAIL 
B26 EPS= 16e**(-13) - 
B26 IF (KNU6EQe9) GO TO 26 
B2P7V6 99 28 J= 1eKNU - 
B26 TOMS NUCIJ) 
S296 CALL UNPACKIIOD) 
AAs CALL FTRAN(I) 
S316 oN 25 T= 16 
332. 25 GiCTeJI= =— VOT) -- 
BIB. 28 CONTINUF : 
4346 PF CALL SHWIFTRGC1 6 3) | 
32%, CALL FTPAN(1) al 
QI. 99 22 t= 1M . 
337. 22 AAtIT)I= YOL) a 
ABS 39 TF (KMU6CEI2O0) GN TO 149 
S346 SO 130 J= 1eKMU 
Aan, TON= MUCJ) ari 


y 
110 | 
| 


TF (TOD eLo eM) GO TC 70 
ZETA- VARLABLE 1S ENTERING 


T8= KINVP ASC TOD) 
1) 34. te to 


> 
a 


YVOONRIVIDRMYVAD 
e) 
4 
us 
iS] 
e 
° 


9 =0 t= 1M 
CIMMASCIDVAEQ 

OATs OTe YVCT) 

FONT NUE 

G2CIIFAS(K DoJ 

CONTINUE 

GO TO 13¢ 


PI- VANIASLE 1S ENTERING 
B= KINKAS(IOD) 


d= 


1 

4 
¥ 
Y 
" 
Cc 
> 
! 

3 
| 
! 

c 
D 


OP rUCvaACsranas 
ac 


<2 
no 
3 


CIPMBAS(T 
NTs OOT+ ¥ 
CONTINUE 
G2ZCITWAS(K) » 
CONT TNUF 
> CONTINUE 
90 150 I= 1M 
Y(t) = 9-20 
TF (KTINS 4S 0NM) eNi0 0) YCKIABAS(RM) D= Le 
CALL ATRAN 
CALL SHIFTR( 3.2) 
HO LAC TS bem 
TF (TOP ASC TINE CT) BACTICEASCI)D = XCT) 
YD 17C J= MPL »NM 
ITF (INFAS(JI4EQ60) GO TC 170 
CALL UNPACK( J) 
ONT= 920 


BRAN 18 2 UN ie sam ye A Me 


DAVIDAANN 


ANANANA™ 


165 
170 


O09 1695 T= 19M 


(LOBAS(T)eEQ20) GO TC 1€5 
2OT= NOT# XCLI*VCT) 
7 ONT I NUE 
BBC IDBAS(J))= DOT 
CONTINUE 
RETURN 


SUBROUTINE ENOPNT( JS ePCl1 oI Se ZFLAGe NET) 
IMPLICIT RFAL*R (4-H»O-Z2) 

RFEAL*A MIN eMIN2 

REAL *4 €C( 809) 


INTEGER PO +SPOLePOZ2 eSeReSSeARe ZFLAGeRS oP 

INTEGE2*2 JH4(350) sO1GMA (952) eK INBAS (1302) » LOBAS( 1302) 
INTEGER &> TSTYPE sUAsLE oI As IE ePUNsLC(20) 1 C1600) 
COMMON/NEWTZ HELO e114 eo XE10 2e ZE10 2s ACC 3e 1029 BLAME 10) 
COMMON/LO1/P101392)eXxK( 13202) 

ZOMMON/ ELTST/ BFL (10 )s BE2C1C De ELC10 2, E2€10) Ce IC oLC 
TOMMONSBLCE ST2/010610 010) 0209910) o0F1 (9010) oF2(9010) 
COMMON/SL NCONS G61 (350010 )2G62(400410)-6Al 350)+88( 400) 
COMMON/INOXLZ NUHC1LO) »MUH( 10) eNUC10)eMU(10) 
CUMMINS IND X27 JH e DI GMA. KINBAS, IDBAS 

COMMON/SCAL/. BT «Ns Jue MELAGs THP 10 P ®D.MP Os KEUNS KJAC 
COMMON/DIMZ THeNo™eKMUe KNU o MPI ONM,ITATL 

COMMONS INTS TPS (30) KIETCKOUNTSISING eKEND 
COMMON, TOLF 24 TOLFZ+TOLEL eT IL CVeTHETAs STPMX, STPRO 
DIMENSION UCB) eVICLOde VEZC10) FC 10) OOT( 4) sQHS110) «UL (10610) 
SYIMENSTON G(10) 


THIS SUBROUTINE WILL FINO THE OPPOSITE ENOPOINT OF 
G¥*(-1)(0) wiTH RESPFCT TO TRE CELL DEFINED BY THE CURR- 
ENT SASISs WHERE 


YELTA 
TWIEL=S VScd*DELTA 


MFLAG=9 WHON WE ARE SOLVING FOR THE QUADRATIC APPROXIMATION 

MFLAG=1 WHEN WE ARE SOLVING ONE OF T HYPERPLANE SUBPROGLEMS « 

MFLAG=2? wHEN ONE OF TrFE LINEAR CCNSTRAINTS IS BINDINGe 

MFLAG=3 WHEN THE FQUILIBAIUM POINT APPEARS YO BE IN THE CURRENT 
CFile 


GOK eTHETA)= THETARFE(X) = CI-THETA) *X 

F(X) = THE SUDGET SURPLUSES DETERMINED BY THE 
FIN AND SLC INS SUBROUTINES>s 

THE TA= THE HOMOTCFY FARAWETFR,. 


eOl 


NEWTON'S METHOC wItt PE IMPLEMENTED AS A TAIL RUUTINE ON Fe 


MFL AG=0 


] 
q 
a | 
] 


TTPY= 0 

THP1l= [He 1 

© (THPL) = 029 | 
KMUt= KMUt 1 . | 
IcK= THPY | 


INTTIALIZE THE INDEPENDENT VARIAGLES IN Xe f 


(KYU 6809) GO TC 615 
YO 8619 T= LeKMU 
$10 x(Td= PL(Mutt)) 


112 


+ ne 


S 
° 
. 

2) 
Cc 
7 


1 
SUM# vV2(1TK)*GI(JE.T) 
105 CONTINUE 


TE (SUM el T2000) GO TO 140 
GO FO 135 
117 sSs 108450 5S) 
T= (ste KECO) GI TO Ile 
IF (V2CT90eGT69) GO TO 135 
Gb TO 149 
112 90 115 T= teKMy 
SUM= SUM# V2(1)*G2( JB!) 
115 CONTINUE 


9 
. 


b. 960. AS I? pat tet ted | GO TO 62S 
le 090 6?C T= 1+eKNU 
: SOP. KM= KMU¢ T 
f W626 S2C K(KM)= XKCNUCT)) 
Wea. XC THPL) = THETA 
WeRe TTFR= C 
Whe 425 39SGCTs ¢ 
QA7. S (Tek= 1LTT9e¢ | 
WOFee CALL GTN JeGeXeF) 
G5Ae CALL ME RIVG(O eX sGeF) 
970-6 7 1sGcT= 1SGCTe 3% 
Wie 1= (315GC%.GTeIH) GO TO 300 
9726 29 10 T= belt 
O73. PHS CLT) = -HOTe TICK) 
Q74. 19 CONTINUE 
QTE. TE (1CKeEQeIHPL) GO TNO 35 
We 99 24 J= 1CXelr 
o77. JOL= J+ 1 
GTA» OO 23 T= elt 
B76 23 HCT.JI= HC TesPL) 
GRC. a 24 CONTINUF 
«ele ¢ 
Ze c Sotuve THE LINEAR SYSTEM TC FING TEE TANGENT TO 
IAN. ¢ THE CUSVE DEFINED SY G**- 100) © 
Rhee € 
2A. 35 CALL JE COUMPCLH.H4eUL? 
GRA, 37 LF CIE INGeNE*1) SC TI 28 
SPT. IF (YCKX eFNe 1) GI TO 35 
WARS “xe FCxX- 1 
aad, Gy tH ? 
7OC, sa Icx= IM 
791-6 or to 
IWWe SA CALL SOL VE CI 4 eUL eRHSe v2) 
99 34 T= 16TH 
IM= THOE1- 1 
2956 IF (iMelL TeICKX) GO TO 41 
IDAs Tei= (Me 1 
IT. ve(TePr1r= v2cImM) 
aan. Oy CONTINUE 
QV365 41 V2CICK)= 19 
We TF CITE ReGTel) GI TO 45 
OOle SUM= 760 
TI32e IF (PUef Qe-1) GO TO 110 
99% 6 J@= KINBAS( JS) 
NCOs TF (JP eNF&2O) GU TO 104 
3086 TF (V2C1996GTeO) GO TO 135 
OOF « GO TO 14¢ 
IC? eo 104 99 108 L= 1eKNU 
D90F 6 Ik= KMU+ 
9 
3 
9 
9 
0 


ee ee eee en ee 


VvIVWVew 
a ee et ee 


A WP at ere > tp ey 


ORPAF UWE eV 
eeereetrteere 
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192¢. e€ 

1921.6 ¢ INIT IS THE CURVE INODEX$ IT TELLS US WHICH OLRECTION 
1922.6 € THE TANGENT MUST PIINT TO KFEP GOING IN THE SAME OIRECTIONe 
1023. Cc 

1024, TF (SUMeLTe O20) GQ Tu 140 

1025-6 135 (NIT= KDET 

L132? he 6Q TO 45 

1327. 1490 INIT= -KDET 

13278656 435 CALt NORM(V20S1 + fHP1 ) 

102%. Sis (KOC Te INITIZS1 

193. 7 47 T= Le THPl 

19316 ATC(2e¢T)= v2(1)*S1 

10326 ACC (3eTd= X(T) 

1933. 47 CONTINUF 

1934. WRITE (462889) 

193F. DN 48 J= lelHPl 

193. WP ITE (52855) Fl ID 6 CACCCI oJ) oe T=2,3) 
1037. 4&2 CONTINUE 

103%. Cc 

LD3Ge S FINC THe FIRST BUUNDARY THAT THE TANGENT. GC ALPHA) 
104°, c 4ISe wr ARE TEYING TO FINE 

1041. ¢ THE SMALLEST ALPHA SUCH THAT 

1942, ¢ Glee *Q( ALPHA) *# GE(T) = 0 

194%, ¢ 

1044, IFLAG= © 

L134, IF (KNU«EQe0) GU TC 249 

1946, 99 25° T= Lem 

1047, lv= JHCOT) 

1048. TF (IT VeFQeM eORe [Veleelh) GO TU 2E0 

1049. 09 240 Lie 2,2 

195Ce U(LLI= O40 

1051. DN 230 K= 1¢eKNU 

1982. K5J= Ke KMY 

198%, VCLLI= UCLLIF GLCL eK) FACCILLEL KS) 
1954. 23? CONT INUE 

1958. zac CONTINUE 

1756. 241 UC3)= UC3+ BAI) 

1957.6 TF (DAAS(UC2Z)) eLTe-TOLFZ) GU TO 250 

105. ALPHA= -U(3)7U(2) 

19596 IF (ALPHA eL Te =TILFZ) GO TG 250 

19606 TF (ALPHA e¢GTs TOLPC) CO TO 244 

L741. TOO= 960 

1%62. 90 242 K= LsKNU 

19636 KM= Ke KMY 

1966.6 242 TOO= TOD Gill eK) *ACC(2]eKM) 
1656 IF (TOO eGTe -TOLFZ) GO TO 250 

1964. 2744 IF CTFLAGe0EQel1) GO TO 243 

10¢7. TFLAG =f 

1346". Js= 7 

1364 MIN= ALPHA 

1070. MPO= 1 

19716 40 Til 269 

1977. 243 T= CAL PHASGS eMIN) GO T) 250 

1072.6 Js= I 

TV7TA.® MIN= ALPHA 

107%. MPO= { 

1074, ae PAN TINUE o 
1977. 2749 TF (KMU CFA) GO TO 649 

LO7A. 3 2969 12 Fer 


TF COT GMACT) e&Qde 4) GC TO 259 


i 
114 } 


DQ 256 tle 
U(LLI= 920 
90 254 K= 16KMU 
U(LLI= UCLLIF GACT eKIFACCI(LL SK) 
CONT INUE 
CONTINUE 
UC2)= UC3)I+ FST) 
IF (DABS (UCF)IeLTeTCLEZ) GO TO 259 
ALPHA= -UC3) 7402) 
TF (ALPHA eL Te -TULFZ) GU TO Z2E9 
I= (ALPHACGTe TOLKBD) GC TO 258 
TID= 20 
79 282 K= YeKMU 
YUD= TOO ACCO[]cK)I*G2CT + K) 
CUNT INUE 
IF €TUO 2GTe -TCLFZ) GC TO 259 
CTFLAGeEQe1) GI TO 253 


W946 


Tr CAL MHA SG e VIN) ¢c T3 259 
JJ= 1 
MIN= ALPHA 
MPO= -1 
CONTINUE 


Oe re ge rt ot ne ee Oe 


2eO2%e INFLCEQel) GC TC 251 
6°59 


=USMAT (Cixe® we ' CISTANCE "5F1Se0e* ALONG THE APPROXe te Ze 
¢ THE CURVE HIT * TYPE CCASTRAINT, NUMGER®.I3) 


enh WO BO ee OO pee OS peo bed OB bes * 


SET THE YNITIAL K FQUAL TC QCALFHA) AS AN APPROX— 
IM4a7TTON TO THE ENDPIINT SOLUTECNe 


£eA35) MIN» MPOe JJ 
“NelLTe STOMX) GD TO 255 
MINX STPRE 


= ) 
(H+635) MIN» MED. JJ 


ACCA eT) + ACC(341) 
3 XCT)*ACC(2AeT) 
CONTINUE 
WRITE (64292) MELAG 
FORMATS 41Xe* MFLAG= %e14) 
T©= (MPD .FQe -1) GC TO 2S3 
TE CJHC II) 0e8& Qe NM) GO TC 295 


Oe ee pe te ee ee te 


9 
2 
fa 
9 
5 
1 
1 
1¢ 
1 
t 
1 
1 
1 
1 
t 
1 
1 
1 
1 
\ 
! 
1 
} 
1 
| 
1 
1 
! 
1 
1 
1 
1 
." 
1 
1 
1 
$4 
1 
1 5 
1: 
1 
1 


Se 


F CALL GTNCMFLAG oGe 
CALL NOUSMOGeS1e TH 


aCe WRITE (e935) Te (XCIT) eI Sle lHPl) 
) ale WRIT= (A 4935) Sle (GCS) sJaLsIhPL) 
a2. CALI DEF IVGIYMFLAG eX eGoF ) 
a36 If (Stel TeTILCV) GC TI 280 


aa, I1v= [T+ 1 

NET= NET# 1 

IF (ITetfFe19) GO TC 255 
MFLAG= 1 
MIN= MIN/?2, 


Pee 
NO 
eee 


iD 
e 


4 

49-6 TF (MIN ¢GTe TCLAYD) GC TO 285 
RE» sT2P 

Sle 266 IF (NET eL Te 60) GC TI 268 

S26 269 wRIT= (6+267) 

S32. ?67 FOOMAT (¢ ENDPCINT H4S FAILED --- TOO MANY ITERATIGNSe *) 
S4. sT7P 

$5) 2768 CALL DFCOMPCTHP 1! sHeUL) 

S66 TF (YTSINGe fel) GC TO 300 

S7e TALL SOLVECLHP LeU eGoZ) 

5e 00 279 T= lel HP 

RO. 2790 XCI)D= XCT)=— 201) 

65 GO ™. 265 

€l. 280 ITF (XCTHPIL) eGTe 1460091) GC TO 295 
62. TF (MPL AGCEQ2 2) GO TY 800 


GO 79 625 


an 
Pa 
ee 


IF CETATL oF 205) 4) TO BO 

TALL OSTNTCALPHe 51) 

IF (ALPHeGTe TOLZV) GC TO 70 

WRITE (49945) | 
50 T2 3¢¢ | 


7° WRITS(62955) ALPH | 
#9 od 9 2 tel 
X(L)= XCD) = ALFH*Z(T) 
95 CONTINUE 
ar 


we 


CONTINUE 
52 72 306 
WE CALL CONCHK(GMING KGMIN» MP2) 


et ee ee ee ed ee ee ee ted ed el ed el ed al el al 
ee ee ne ed deen ed ed et ed ed ee dd el 
> 
e 


- 
ES. 9 NEWTON ES METHOD FCR THE PROBLEM CF FINDING THE INTERSECTION 
Fe Cc CF THe CURVE WITH THE CONSTRAINT DEFINED 38Y MPO AND JJe 
a, e 
295 MFLAG= 3 
69 90 19¢ K= 1.221 
We *¥M1=e K= 1 
Tle CALL FYN( FX) 
72 CALL NIU2MOS 6 SlelIH) 
736 WRITE (6e93A) KML oe (XC Te l=zlelIb) 
74.6 WRITE (66939) Slo CFC(L) oL=1leTH) 
7Se TF (SLeLTCTWEFEZ) GI TO 800 
TF. CALL DERIV 
77. NET= NET+ 1 
The IF (THeGT.et) GO TO €§ 
7%. IF (DABS(H(1+1)d) el Te TCLFZ) GO TO SB 
BP. ZCLI= FCLIZAC 161) 
Ple GO F > 65 
P2e Ss CALL D=7 IMP (1k erty UL) 
F&3. TF CISINGeFQet) GD TO $C 
R4, ba) WRITE (6,949) | 
CE. sToOPe : 
B5e 450 CALL SOLVE CTHeULs FeZ) 
7. 65 ALPH= le 
a 
R 
9 
9 
S 
9 
co) 
9 
y 
a 
g 
9 


UD PNe ArI—MPFITVF NP wD 


erteeterteere 
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Ft2 SU STO PAR ee stm ho 
(i el age a APL Si SRL aa a Rn 


i—_— — oa 


THE YA= X(LTHYP1) 
TF (KMU6CEGe9) GO TC 146 
00 145 [T= 1.«kMU 
KM= MUCTI) 
145 OTCKM)= X(T) 
146 IF (KNUCEQNSS) GL TOC 1438 
99 '47 T= KMU1L eH 
IK= I~ KMU 
KM= NUCIK) 


120C. 

1201. 

12026 

12036 

1204. 

1205. 

1206-6 

1207. 

1208. 

12CGe 147 XXCKM)= XT) 

121%e 1428 TF (GMIN eLTe =TILAN) GY Te 170 

L2lle IF (MFLAGeLE.2) GC TO 133 

1212. 7ZFLAG= 2 

LPLd.e RETU?N 

1214.6 1S3 IF O(MFLAG.EQ.1) 30 F? 159 

1215e FOS FE €4PO erGe ~1} GO TT? 181 

L216. TS= JHC JJ) 

1217.6 PNI= 1 

1271, RE TUN 

1215.6 181 §S= OTrGEMAt I) 

L2?%. POL= ~3 

$226. RE TURN 

12726 139 MFLAS= ¢ 

17236 GO E125 

1224, Cc 

12756 € IF THE CURRENT VALUF OF X VIOGLATES THE CONSTRAINT DEFINED 
1224.6 1S SY KGYIN AND “OP, FINC A COOD STARTING POINT FOR NEWTON'S 
1227. e ME THU 2Y FINO'TNG TH: POUINT CN THE LINE SEGMENT 
12°. « ACCOLese ee) + ALPHA*(X - ACCOLleed de O<= ALPHA<= 1 
12256 a THAT 347TISFIES THE VITDLATEC CUNSTFAINT EXACTLY «© 
123% © ¢ 

eal 179 JJ= KGMIN 

1237 MPD= “WP 2] 

¥23 30 MFLAG= 2 

L234, WRITE (65174) YEDeJI+GMIN 

1 72% 's *74& FORMATOCZ,* THE QUACRATIC PICKED THE WRONG CONSTRAINTe's7Z 
112A. 1° THE MOST INFEASIBLE CCASTRAINT 15 TYPE *oI3¢ 
T2376 Z2* NUMSER "elAs/s* WITH A VALUE CE * oF 1407) 
173-6 99 £71 = ieIHPl 

123%. 171 VICTI= XCTI— ACCO3e1) 

124C. OD 3175 [= 164 

1241. Lets NOT(T)= 920 

1242. IF (“PDeFQe-1) GO TH 189 

12423. 00 172 T= 1¢KNU 

1244, Ik= [4 KMU 

124°. NOT(L I= TITCLI+ GICIS eo T)RACCCZeTK) 

124Ff 6 IUT(2AV= CITOCAI4 SICISeoLI*VICIK) 

1247 6 iv? CONTINUE 

124A STS= (-8aAC II) -9NT(1IIZ CITC] 

12476 WRITS CAZ1RBS) GICISol?) FAC UI) eOCT(C1) e9CT(2)6STR 
12S5Ce 99 173 1= Let 

12516 “73 XCT)= ACCOIST)+ STRRVICT) 

1257.6 690 TY) 26a 

125.76 136 O00 13°72 [T= 1K 

1284.6 MTEC LP= ONTCVI+ G2CSIeTI*FACC(36T) 

32°56 COTC2? I= NITC2V4+ GPCIFeTIFVICT) 

L2ASF. 1°3 CONTINUE 

1257. STP= €-9805I)-DITCLIISOCT(CE) 

1 25e, WRITS CO 6185) GPC S5o1) EEC II)» TOT(1)-90T(2)¢STR 
tens YAS FORMAT (Xe ® G20 BSe BDOTLs DOTA 6 STR 2% -SFI3e6) 
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l 
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1 
3ZeT+ STx*VICT) 


CURFENT QUADRATIC APPRUXIMATION TO THE CURVE IS?* 
x) A2C*UDe A3*) 
17e7) 
ERATION *e14e* X= %, GF1IO06S) 
RMCS(XD)= SsF ll eBs® G= *,3F11.6) 
SP aRKAL GIR ITMM BONBED wITH SINGULAR JACOBIAN***** ) 
ITS RATION 2 %elSe® X= "5 9F1045) 
*"NORMCF (XD = "ePbleEs* F= *eIF100 6) 
FORMAT *NO DESCENT PCSSIBLE -= ALPHA WAS ZERO. ') 
EF OQRMAT ( STEP STZE= ALPRKA= *»FG67) 
FORMAT ¢ "NEWTON'S MSTFOD FAILED TO CONVERGE ce? ) 
gsT qe 
FNS 
Qu40S tFINOS THE SMALLEST NUNNEGATIVE &RJOT OF 
UCT *ALOHA "2 + UCE)*FALYPHA +€ UC3) = Oe IF THERE 
ES ONF e L© NOT, TVAG I35 SET 2 Be 


ADBPIPAROD 


> 
PN: 
ee 


NwwNA 
TRE iA 
eere 


BANN 
aA at N 


~ 
e. 


SUPROUTINE QUADS (Us IMAG sAL PHA, BETA) 
IMOLICIT Re aL*® (A-HeD—2) 

INTEGTS BP) ePOL eH IZ] eSeP® eo SSeRReZFLAG sRS oP 

INTE GEE? INC 350) eH IGYA(952 eK INRAS( 1202) - IODBAS( 1302) 
COMMONS TOLE RZ TOLF Ze TOLIC e TCLCVeTHET Oe STPMX,STPRO 
ZTOIMMONSOTAZS THeNeMeKMJ eo KNU oMP IL oNMolI TAIL 

SIME NSTON U(3) 

IMAG= ¢ 

IF (o4eNs(y TOLFZ) G3 TO 10 

IF (OA32S¢U TOLF2) GC TO 8 


VNU VYVIVVN URN NV NINN ANN VIR NNy Viviv 


Qe TYSN 
UC2YKE!] = 42%U(1)*UC2) 
CET} 209 306 40 


TREWVFDODV ADP 
eeereerteeeetetee eteee 


ODNOIIOVO MDT p CH LO 


1 
1 
1 
1 
1 
1 
1 
b 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
bf 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 
1 


NivVinivw 


THE QUADRATIC NIES NIT INTERSECT THIS CCNSTRAINT. 


IMAG= } 
RETUUN 


THETE 15 A DOVUBL = FOUT. 
ALPHA = UL2)/(UC1IF20) 


QETA= ALPHA 
ae yuS NJ 


ao. 
ss 
ve 
fe 


Two =) AL BIOS 


NIGS= MQRT CRT) 
ALPHA= (-U(2)=D 
BETAS (-U02)401 
TF CALPHA oLEo 
SAVS = ALPHA 
ALPHA = OTA 

Bt Tas SAVE 
SSTUEN 


Yo MW we 
CeryarnFeure>d 


ee er er 


132C. ENO 
L32Qte C 
13226 SUARNUTINI DSENT CAL PHA e-ENGRM) 
L324 6 IMPLICIT REAL*3 (A-HeI-Z) 
Lt?a. COMMOANSZNGE WTS HOLDeL DL eXC1LC Ve ZE1C 26 ACCO3410)¢ RL AM( 10) 
13256 T IMM INSOUM ZS THe No Me KM eo KNU MPL OAM OT TAIL 
L326 YEMENCTUN VOICE )2F (19) 
L327 « $ 3A 6 3s. fai 
132" » > YCTI= XCTI— ALPHA*Z( TI) 
1429 CALL FIN(F,.Y) 
1430.6 CALL NOEM (F5S2eTH) 
Et ie TF (67 el Te PNORM) RETURN 
1732. AL PHA= ALPHA/?,. 
LIF TE CAL SHG «GT « 2 3OO9)) a TD 2 
L334. SEU 
13356 eNO 
1%3ce tS 
1337.6 SURFRUTINE GTINC IRL eQsUeF) 
34386 IMPLICYTY S£AL #3 (AH, 3-2) 
13296 INTOGTF 2D ePO1 sO 92 0Se xe SS5e RR oe ZFLAG PRS oP 
234°. INTEC ~F > JHC 250) eo IGHMA( GSES?) eKIR-AS( 1 392) e LO5AS(1302) 
1341. TIMMONS NiWwTS HOLD C eT Le ACI De ZC1C) +e ACK(C 2+ 10) oe BLAM( 10) 
1 $62, COMMER SOEWS Tet eNne eK e KNUe MPL eAMSITAIL 
1347.6 COMMONSE CALS BT NBeJIS eo MEL AG oe LTHA1 oP oP eg AOD oe KF UNe KJAC 
1344, COMMON SENC ONS ZGI1 (352015 )2G2(4590 610 )e3A(355)eBR(40C) 
144s, TOMMOUNSINOX2Z JH eDIEGMA eKINBAS Ss IOBAS 
1446, TAMM ONS INO XEA NUACLO) + MIRO 1096 NUC1C) eMU010) 
1347. DIMENSION UC19) -6-9(19) oF (10) 
134", QC THE y= OC 
154%, CALL FINCF oJ) 
L350e THE TA= UCTHOL) 
L351. TF (KNU6=O060) GD TC 21 
13526 90 27 b> 2e6N0 
t 3536 Ik= T+ «¥U 
1384. QO5= ©CIK)+ UClK) 
13556 Q(t )= THETA * 9S = UCIK) 
13566 et AONT INGE 
¥ 2875 7. fF €£8Ue-QeGF GO FC 24 
1L%55 6 DV SF £= 14KMY 
13596 Tos MUTE) 
1 2450 o Ide = KINBAS( IP) 
13616 SUM = 909 
17626 T= (KNU e&Ge O) GJ TC 7 
13436 > 2 K= LeKNU 
13646 "K= Kt KMU 
by ee 3 SUM= SUM# GI(ITO3eK) *® UCT K) 
136éFe ve ALAIMCT = SUM + EAC TER) 
1367 6 TFs FLAMCT) 
L3EF° Qs FUT)* OFF 
135°. Q¢t)>+ THETA*®3S - OTF 
he Sab 23 CONT (INUF 
L374. 24 TF CIEL SFG. GO} KETURN 
1375.6 KF UNS KEUNe 1 
13766 sus ier 
1377.6 Be OT" Ue Oes) Goarrke sy 
3.3 ee wy 2S T= Te Peery 
137°. SUN = SUM + ACC(AsT) * UCT) 
1780. 2 ZTONTINUF 
13Ale QCTHES = SUM = AT 
17982, 9e TURN 
= 
ei 
ie 119 


2 
ud 
e 


ITF (MP) .FQe21) GY TI 5C 
IF (KNUSFQ29) GO TC 45 
09 40 T= 1LeKNYU 
Ix= Tt KMU 
SUM= SUM+ GI(JJeoTD FLCI KD 
“ONT TNUF 
QCTHP 1) = SUM+ BAC JJ) 
wE TURN 
TE CHEMO’) GO TD 49 
AN SD T= 1e¢KMU 
SUM= SUM+ G2CISeT) UCT) 
CONT TNUF 
aC The 1)= SUM+ SBC JY) 
RETURN 
=N9D 
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SURTIVTING OF RIVG CALCULATES THE HESSTAN OF THE 
RTT OACTT CS, FUNCTIONe PLUG SCME OGTHFR GRAUIENT 
SPECIFIES AY IF s« 
Teos0 NO OTHER GRACTERT 
TFL=t 4 SPHERE DITcRMINES THE GRADIENT. 
TFs? ONL OF THE LINEAR CCNSTRAINTS IS 
THE GRADTENT. 


SUBEMUTINE DERI VGC IE LeU eGoF) 
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